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CHAPTER 1 


INTRODUCTION 


The displacement formulation of the finite element method is the most general and 
most widely used technique for structural analysis of airplane configurations. Modem 
structural synthesis techniques based on the finite element method have reached a certain 
maturity in recent years, and large airplane structures can now be optimized with respect 
to sizing type design variables for many load cases subject to a rich variety of constraints 
including stress, buckling, frequency, stiffness and aeroelastic constraints (Refs. 1-3). 
These structural synthesis capabilities use gradient based nonlinear programming tech- 
niques to search for improved designs. For these techniques to be practical a major 
improvement was required in computational cost of finite element analyses (needed 
repeatedly in the optimization process). Thus, associated with the progress in structural 
optimization, a new perspective of structural analysis has emeiged, namely, structural 
analysis specialized for design optimization application, or what is known as “design ori- 
ented structural analysis” (Ref. 4). This discipline includes approximation concepts and 
methods for obtaining behavior sensitivity information (Ref. 1), all needed to make the 
optimization of large structural systems (modeled by thousands of degrees of freedom and 
thousands of design variables) practical and cost effective. 

In the airplane conceptual and preliminary design stages configuration shape optimiza- 
tion is essential. Wings should be allowed to change in planform and airfoil shape. Fuse- 
lage structures should be allowed to be shaped simultaneously, and the position of wings 
and control surfaces should be determined as part of the optimization process. While a 
substantial amount of work in the context of structural optimization has been devoted to 
structural shape synthesis of solid and machine parts, very little has been done to date in 
the area of aiiplane structures. Moreover, even with the availability of computer graphics 
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and computer aided design tools, the preparation of a new finite element model for a new 
configuration is still too time consuming. It is estimated in Ref. 5 that it would take about 
12 months to complete a single structural, loads and aeroelastic design cycle for a high 
speed civil transport. A major part of this effort is dedicated to the generation and updates 
of the finite element model. 

This thesis focuses on techniques for modeling airplane wings for the conceptual and 
preliminary design stages using finite elements. The emphasis is on shape optimization. 
An automatic mesh generator is developed to efficiently handle planform and airfoil shape 
variations. Simple bar and triangular membrane elements are used to represent spar / rib 
caps as well as skins and internal webs. Analytic deformation, stress and natural frequency 
behavior sensitivities are obtained with respect to shape design variables in addition to the 
sizing type design variables. Extensive numerical tests of the resulting modeling technique 
are conducted to evaluate its accuracy and economy. The new technique combines advan- 
tages of equivalent plate wing modeling (Ref. 6) (ease of model generation and shape sen- 
sitivity calculations) with those of finite element models which are general and can handle 
local effects and structural discontinuities in real wing structures. 

The outline of this work is as follows: in Chapter 2 the two finite element modeling 
approaches are discussed. In Chapter 3 wing behavior sensitivities with respect to both 
shape and sizing type design variables are derived. Chapter 4 will focus on aspects of 
automatic mesh generation while Chapter 5 will deal with issues of finite element model- 
ing implementation. In Chapter 6 the three wing models to be analyzed are introduced and 
described. Chapter 7 details all results pertaining to wing displacements, stresses and nat- 
ural frequencies while Chapter 8 concludes with sensitivity and computational cost 
results. Detailed mathematical derivations are given in the appendices. 



CHAPTER 2 


MODELING CONSIDERATIONS 


2.1 Introduction 

Two modeling approaches for built up wing structures are described in this chapter. 
Both are based on truss (rod) elements for spar and rib caps. Membrane (plane stress) ele- 
ments are used for cover skins and spar/rib webs. The motivation for using these simple 
models is not only in their simplicity and speed of computation, but mainly because it is 
possible to obtain closed form explicit analytic sensitivity of their stiffness and mass 
matrices with respect to shape design variables. In the first approach linear rod elements 
and constant strain triangular membranes (CST’s) are used. In the second approach linear 
rod elements and linear strain triangular membranes (LST’s) are used. Discussion of these 
two approaches and guidelines to follow are included in this chapter. In both cases there 
are no rotational degrees of freedom in the model. 

2.2 CST modeling 

The simplest of the two techniques is the one employing the three-noded CST mem- 
brane element with a linear rod element. The CST is used to represent all wing cover skin 
panels and rib and spar shear webs. The rod element is used to model all rib and spar cap 
areas. These are low order elements. Stresses in these elements are constant throughout 
their interior and for convergence a large number of elements may be needed. 

A finite element capability, then, must include grid refinements that are quick and easy 
to perform, and a study of modeling accuracy to establish modeling guidelines as to the 
degree of grid refinement required. The possibilities to be investigated include refinement 
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in the spanwise direction, refinement in the chordwise direction or a combination of the 
two. For grid refinement in one direction only, more nodes are created along the spars for 
refining spanwise, while more nodes are added along the ribs for chordwise refining (Fig- 
ure 2. 1 ). As one can see, all newly created nodes still lie on a rib or spar, and thus they are 
supported by the internal structure of the wing. For a combination of the two, grid refine- 
ment introduces a “floating node,” or a node that has no vertical support (Figure 2.2). As a 
result, the stiffness matrix becomes singular, and a special procedure must be used to elim- 
inate this singularity. One way of overcoming this difficulty (Ref. 7) is by linking the dis- 
placement at a floating node via multi-point constraints to the displacements of it’s 
neighboring nodes. Since the equations of constraint depend on wing geometry, though, 
analytic differentiation of stiffness and mass terms with respect to shape becomes quite 
complicated. 

Our solution is to add either “dummy” rib or “dummy” spar elements whose thickness 
is substantially lower than the real ribs or spars (say, 1% thick) so as to not influence the 
stiffness or mass of the wing but provide support for the floating nodes. The advantage 
here is that all nodes and elements (whether real or dummy) are treated in the same way in 
the course of analytic differentiation and no special treatment has to de devised for the 
floating nodes. It must be remembered, however, that floating nodes can not have any ver- 
tical loads applied to them. Thus when aerodynamic loads are distributed over the wing 
they can only be applied to nodes supported by the actual internal structure of the wing. 

Using CST webs for the spars and ribs creates another problem. Since only one row of 
CST elements is used in the depth direction of the wing due to a wing’s small depth/chord 
ratio, this leads to finite element models that are too stiff (Ref. 7). This comes as no sur- 
prise since the constant stresses in a CST cannot capture the linear distribution of stresses 
in a typical beam web. To correct for this the CST web membrane elements are modified 
to only carry shear stresses by using just the shear stiffness portion of a CST’s stiffness 
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Figure 2. 1 - CST spanwise or chordwise refinement 
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Figure 2.2 - CST spanwise and chordwise refinement 
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matrix. They act as pure shear webs, and vertical rod spacers connecting the upper wing- 
skin to the lower wingskin replace the normal stiffness of the web elements in the trans- 
verse direction to keep upper and lower skins separated. 

2.3 LST modeling 

Using LST elements in place of the CST elements leads to better convergence of finite 
element results because of the higher order of the LST. The problems with floating nodes, 
however, are still present. The LST is a 6-noded element whose three additional nodes are 
located along the midpoint of its sides. Because of these midside nodes, floating nodes 
now appear not only in the wing skin planes, but also in the rib and spar planes (Figure 
2.3). Thus, in the spirit of our approach to CST modeling, a combination of two of the fol- 
lowing techniques is necessary to provide support for these nodes and eliminate singular- 
ity: dummy ribs, dummy spars, and/or dummy layers. The dummy layers are added to 
support the mid-depth nodes of the spar and rib webs. They are similar to the other dummy 
elements in that their thickness is vei^ low (1% of the actual wing skin thicknesses). 

Since the LST’s lead to better convergence of the finite element solution, the most 
basic mesh possible (the one defined by the location of real spars and ribs in the wing) is 
usually quite accurate. The stress output for a LST element consists of a pair of normal 
stresses o , a and a shear stress a at the comer nodes where each stress varies lin- 

xx yy xy 

early across the element’s interior. 

Due to the higher order of the LST, shear stresses through the wing thickness are better 
represented. Thus, no pure shear LST is necessary when using LST models. 
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12 elements total 
4 'dummy* elements 


Figure 2.3 - LST dummy element selection 
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2.4 Wing lumped mass modeling 

For natural frequency calculations, a lumped mass matrix modeling technique is used. 
The mass of each finite element is distributed evenly to it’s nodes and then merged to the 
global lumped mass matrix whose structure is strictly diagonal. Since floating nodes can 
not support any load or force, inaccuracies in the calculation of higher natural frequencies 
and mode shapes will arise. It will be seen in Chapter 7 that natural frequency accuracy is 
a direct function of dummy element thickness and guidelines for the selection of this 
thickness will be provided. 

2.5 Finite element derivations 

For complete details of the finite elements used and their respective stiffness, stress 
and mass matrices, consult Appendix A. 



CHAPTER 3 


BEHAVIOR SENSITIVITIES 


3.1 Introduction 

Accurate and computationally efficient derivatives of behavior functions (such as dis- 
placements, stresses or natural frequencies) with respect to design variables are important 
in the context of gradient based optimization not only for calculation of the gradients 
themselves but also as a basis for constructing constraint and objective function approxi- 
mations (Refs. 1, 8 and 9). When structural shape optimization is involved, it is difficult to 
obtain these sensitivities in a closed, explicit analytic form (without any numerical inte- 
gration, as is usually used for evaluating mass and stiffness terms of general elements). 
One popular way for obtaining structural behavior sensitivities is by finite differences 
(Ref. 1). This technique, however, can be time consuming when the computational cost of 
a single analysis is high. In addition, and especially in the case of shape variations, finite 
difference derivatives are sensitive to the step size used, and can lead to erroneous results 
(Ref. 1). 

In the present finite element modeling capability developed, simple finite elements 
such as truss rod and plane stress CST’s and LST’s are used not only because of computa- 
tional efficiency in formulating the stiffness and mass matrices, but also because of the 
explicit algebraic nature of these matrices (Refs. 10, 1 1 and Appendix A). This makes it 
possible to obtain behavior sensitivities in an analytic, explicit manner, thus avoiding 
numerical problems associated with finite differences and significantly reducing comput- 
ing time. 

The wing structural design variables are divided into two categories: shape and sizing. 
The wing planform is divided into trapezoids. The shape of each trapezoid is defined by 
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six shape design variables. The variables y L , y R are the left and right spanwise coordinates 
of the trapezoid, while xpL, x^, xpR, xar are the longitudinal locations of its four verti- 
ces (Figure 3. 1). The sizing variables include the cross-sectional area Aj of any rod ele- 
ment T and the thickness tj of any CST or LST membrane element ‘j.’ Based on the 
formulations in Appendix A, analytic expressions for the sensitivity of element stiffness 
and mass matrices can be derived with respect to the location ot an element’s nodes. This 
is done here in a manner similar to Ref. 12. The position of each element’s nodes can be 
linked to the overall shape of an individual wing trapezoid knowing the rules used for gen- 
erating the mesh for that trapezoid. Chain rule differentiation is then used for obtaining 
stiffness and mass sensitivities of individual elements with respect to overall wing plan- 
form shape design variables. Details of these derivations can be found in the appendices. 

3.2 Sensitivities with respect to shape variables 


3.2.1 Global displacements 


The linear static structural equation serving as a basis for static analysis is 


[K] {U} = { F} 


(3-1) 


The equation for displacement sensitivity with respect to any design variable in the 
case where external loads do not change is (Ref. 1 ) 


9{f/} 

~W 


, rl -Kam 

= -<H 


{U} 


(3-2) 


where [AT] is the stiffness matrix, { U} is the displacement vector and P is a typical 
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design variable. 

Once displacements and displacement sensitivities are known, it is possible to obtain 
derivatives of stresses within elements. The stiffness matrix [K] is nonlinear in the shape 
design variables. However, explicit expressions for stiffness terms in rod and plane stress 
elements are available (see Appendix A) and can be used for differentiation. 

3.2.2 Stress in the i’th rod element 

As shown in Section A. 1.2, the stress in a truss element depends on the shape design 
variables both explicitly (through a location vector {X}) and implicitly (through an elastic 
deformation vector {Uq}). Therefore 


dcr 5a. d {X} . ^ 8o ( . dif/c),. 

W = a{x} j . ap + a{f/ G }. ap 


where { X } j and { Uq are the location and displacement vectors in global coordinates 
associated with a rod element, respectively. 

3.2.3 Stress in the i’th CST element 

Stress sensitivities for the CST with respect to planform shape design variables are 


obtained by differentiation of the stress equations in Section A.2.2 giving 
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d_ 

5p 



d{u r ) 


a I a] si«i . _ 

[«], [»],-5T' *C'il.-5r0.' ‘'■■■I, 


(3-4) 


where { Uq } j is the vector containing CST element i’s nodal displacements in global coor- 
dinates. The material matrix [D]j does not depend on the shape of the element, therefore 
it’s derivative with respect to planform shape is zero. 


3.2.4 Stress in the i’th LST element 


The equations of Sections A.3.1 and A.3.2 are now differentiated analytically to obtain 
sensitivities of stresses at the three vertices of an LST with respect to shape design vari- 
ables. Chain rule differentiation is used to link variations in element node locations to the 
global planform shape changes of the wing to give 



(3-5) 
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The matrix is a material constitutive matrix and does not depend on the shape of the 
element. {Uq} is the vector of element nodal displacements in global coordinates. 

3.2.5 Natural frequencies 


The governing equation of motion for an undamped structure in free vibration is 

[K-(0 2 M] {<D} = {0} (3-6) 

where X = CO 2 is an eigenvector, { } is it’s respective mode shape and co is a circular nat- 
ural frequency (radians/second). Implicit differentiation of X with respect to any shape 
variable (3 yields 



for eigenvalue and mode shape ‘i\ Since the natural frequency (in Hertz) is given by 



the natural frequency sensitivity after differentiation is 

df. I dX. 

" 4* 


(3-9) 
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3.3 Sensitivities with respect to sizing variables 

In this case the stiffness and mass matrices depend linearly on the design variables. In 
the case of truss elements and CST’s or LST’s, these design variables are cross sectional 
area and membrane thickness, respectively. 

3.3.1 Global Displacements 

With K as a sizing type design variable, the matrix equations for sensitivities of the 
displacement vector in global coordinates are (Ref. 1) 


at u] 

~3¥~ 


= -<&r> 


-i.am 

3k 


{U} 


(3-10) 


Again, it is assumed that external loads do not change with changes in the sizing design 
variables. 


3.3.2 Stress in the i’th rod element 


If the design variable is a rod cross sectional area Aji 


da i 3a. B{U C ). 
BA; = BTUgT. 3 Aj 


(3-11) 
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If the design variable is a membrane thickness tj: 


9a. 9a ( . 9 {U G } . 

u: = 9 { u G ) . dij 


(3-12) 


where { Uq }j for both stress sensitivities is a 6x 1 vector containing nodal displacements in 
global coordinates for rod element i. 


3.3.3 Stress in the Pth CST element 


If the design variable is a rod area Aj: 



9 { U r } . 


(3-13) 


If the design variable is a membrane thickness tj: 




9{t/ r }. 


(3-14) 


where { Uq } j for both stress sensitivities is a 9x 1 vector containing nodal displacements in 


global coordinates for CST element ‘i.’ 
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3.3.4 Stress in the i’tli LST element 


If the design variable is a rod area Aj: 



If the design variable is a membrane thickness tj: 




-O' 


*y 



(3-16) 
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where { Uq } j for both stress sensitivities is a 18x1 vector containing nodal displacements 
in global coordinates for LST element i. 

3.3.5 Natural frequencies 

With k again as a sizing type variable Aj or tj. the eigenvalue sensitivity is 


9X. 

I 

dK 


Of 


ram , 


9k 


-X.- 


9k 


*/ 


<t>f [M] <t>, 


(3-17) 


and the natural frequency sensitivity is 


3/, = I 
3ic 4n 


Derivatives of nodal displacements with respect to shape type design variables are 
obtained from Section 3.2.1 while derivatives of nodal displacements with respect to siz- 
ing type design variables are obtained from Section 3.3.1. Sizing derivatives of the stiff- 
ness and mass terms are straight forward because of the linear dependency (see Appendix 
A). Thus, stress sensitivities with respect to sizing type design variables require sensitivi- 
ties of deformations only. Additionally, all other matrix and vector transformations used to 
move from local to global coordinates and from deformation (displacement) to stresses are 
fixed in this case. 



CHAPTER 4 


AUTOMATIC MESH GENERATION 


4.1 Introduction 

The desire to circumvent the creation of finite element input files by hand and to auto- 
mate model generation for wing shape synthesis makes it necessary to combine a mesh 
generation capability with the finite element analysis and sensitivity techniques (Ref. 13). 
At this stage of the present work this capability is limited to wings with ribs parallel to the 
root rib, spars beginning at the root rib and terminating at the wing tip and a thickness dis- 
tribution symmetric about the wing’s mid-plane (Figure 4.1). This modeling is sufficient 
for the studies conducted in this work. The limitations are minor and can be removed by 
making the mesh generator more general for other wing layouts and also for fuselage 
structures. For the structural wing model the elements used include constant stress rods (to 
model cap areas) and either CST membranes or LST membranes (to model wing skins and 
webs). The mesh generator and finite element capabilities are linked together so that when 
combined with an optimization package, the shape of the wing (in addition to cap areas 
and skin thicknesses) can be optimized. 

4.2 Wing design variables and design rules 

Figure 4.2 shows a sample mesh created by the mesh generator, and will be used to 
define key parameters needed. The structure shown is a single wing trapezoid containing 
five structural spars and six ribs including a root rib. In order to refine the mesh it is possi- 
ble to add “dummy” ribs and spars between structural ribs and spars, as shown. The 
parameters “adrib” and “adspar” define the number of added dummy ribs or dummy spars 
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>- 



Figure 4.2 - Mesh generator data input requirements 
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between adjacent structural ribs or spars. These dummy ribs and spars, whose stiffness is 
negligibly low compared with the actual structure, are needed to support the added “float- 
ing” nodes on cover skin surfaces modeled by triangular membrane elements. This is nec- 
essary since there are only displacements and no rotations associated with each node, and 
since the skin cover elements have no bending stiffness. 

The x-y coordinates of the vertices of the wing trapezoid and span wise and chordwise 
locations of ribs and spars serve as shape design variables for the planform. All dummy 
spars or ribs are assumed to be evenly spaced between real spars or ribs. Wing depth defi- 
nition is also used based on associated shape design variables. Finally, all spar and rib cap 
areas and all wing skin, spar web and rib web membrane thicknesses are used as sizing 
type design variables. It should be mentioned again that at this stage of the present work 
rib generation is limited to ribs that are parallel to the root rib, and spars have all to origi- 
nate at the root chord and end on the tip chord of a trapezoidal section. 

4.3 Planform expansion to three dimensions 

With wing depth specified by the proper shape (depth) design variables, an explicit 
equation for depth distribution as a function of x and y is established over the wingspan. 
The planar mesh described in the previous section is now projected upward and down- 
wards to generate the meshes for the upper and lower cover skins. Realistic thickness and 
camber distributions can be modeled by proper selection of depth shape functions and 
construction of a series in those functions whose coefficients serve as shape design vari- 


ables. 





24 


4.4 Shape variable coordinate linking matrix 


With the 3-D grid complete, the linking of each node’s x-. y- and z- coordinates to the 
six planform shape design variables of a wing trapezoidal section Xp^, XaL’ x FR- x AR- yL 
and y R (Fig. 3. 1) and its depth design variables is straightforward, as detailed in Appendix 
B. Derivatives of each nodal location with respect to each shape design variable can easily 
be obtained. These derivatives, used in the finite element program for shape sensitivity 
analysis, are the same as the coefficients that link each node to the shape variables since 
the linking equations are linear. 

4.5 Finite element placement 

Individual finite elements are placed and connected to the proper nodes according to 
the following rules: spar and rib caps (for the real structural spars and ribs only) are repre- 
sented by rod elements connecting nodes on the upper skin or lower skin along the spar or 
rib lines. Intersections of spar lines and rib lines (including dummy spars and ribs) define 
quadrilateral cells on the upper and lower skins. Each of these cells is divided into two tri- 
angular elements. For the webs of all ribs and spars, each quadrilateral cell (defined by the 
end nodes of the upper and lower rod elements associated with the cell) is divided into two 
triangular elements. 

As discussed earlier, mesh refinement involves the need for dummy ribs or dummy 
spars if a floating node is present. For CST models, dummy ribs are sufficient. For LST 
models, dummy ribs and dummy layers are necessary to support both vertical and horizon- 
tal floating degrees of freedom. The dummy layer (of negligibly thin material) connecting 
the mid-side nodes of LST used in spar and rib webs is covered by triangular elements in a 
manner similar to the cover skins. 
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Following the rules described above, it is possible to generate explicit relations 
between each element, its end nodes and the global shape design variables defining the 
shape of the whole wing. These relations are then used for obtaining analytic sensitivities 
of stiffness and mass sensitivities using chain rule differentiation (as described in the 
appendices). 



CHAPTER 5 


FINITE ELEMENT IMPLEMENTATION ISSUES 


5.1 Introduction 


Implementation issues concerning the finite element modeling technique described in 
Chapters 2-4 are discussed in this chapter. A standard displacement approach is followed 
(Refs. 14, 15). The finite element code of Ref. 14 for three dimensional trusses serves as a 
basis upon which the new capability is developed. Constant strain triangular elements 
(CST’s) and linear strain triangular elements (LST’s) are added to the library of elements. 
A banded matrix solution solver (Ref. 16) is used for static analysis. A QR eigenvalue 
solution technique is used for the natural modes analysis. Analytic sensitivities of stiffness 
and mass matrices are generated and used to obtain sensitivities of displacements, stresses 
and natural frequencies. 

5.2 Global displacement solution 

The governing equation for a static structural system is given by 

[K] {U} = {F} (5-D 

where K is the banded global stiffness matrix, U is the vector of global displacements to 
be solved for and F is the vector of nodal loads. The decomposition technique of Ref. 1 6 is 


used for solution. 
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5.3 Natural frequency solution 

The governing equation for a dynamic structural system undergoing undamped free 
vibration is given by 

[K-vrM] {<}>} = {0} (5-2) 

where K is the global stiffness matrix, 0) is a natural frequency in radians/second, M is the 

lumped global mass matrix and <j> is a mode shape. For a non-trivial solution of natural 
frequencies and mode shapes to exist, the determinant of [K - 0 ) 2 M] must equal zero. 
The method of solution will be to use a QR decomposition algorithm (Ref. 1 8) that solves 
the standard eigenvalue problem 

[A-\/] {\|/} = {0} (5-3) 

where A is a square symmetric matrix, X is an eigenvalue, I is the identity matrix and \y 
is the corresponding eigenvector. The original equation is converted into the standard form 
by using the fact that since M is diagonal and positive definite, it’s square root is easily 
calculated. Thus, pre- and post-multiplying eqn. 5-2 by (JM) 1 gives 

{\|/> = {0} (5-4) 


or 


[A - \[) = {0) 


(5-5) 
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where A = ( JM) K ( JM) is square symmetric and X = to 2 . 

Since \|/ solves the standard eigenvalue problem (eqn. 5-3), to find § which solves the 
original eigenvalue problem (eqn. 5-2) it is necessaiy to use the formula 

<j> = . (5-6) 

5.4 Element stress solution 

All individual finite element stresses are calculated using the previously found global 
displacements. Equations for axial stresses in rod elements are given in Appendix A. 1.2. 
Similarly, for CST elements, stress equations are given in Appendix A.2.2. For LST ele- 
ments stress equations are given in Appendix A.3.2. 

5.4.1 Stress smoothing 

Since stresses throughout a CST element are constant, stress differences can be found 
between neighboring CST’s and an averaging process (Ref. 17) is needed in order to 
obtain a smooth stress distribution over the skin in areas where no discontinuities are 
expected. 

As an option in the present capability, a least squares fitting procedure is used to fit an 
order polynomial for each skin stress (0**, 0 yy , 0 xy ) over each wing skin trapezoid. 
Thus if S(x,y) is any component of the plane stress in the skin, then 

su?) = <?i + <? 2 * +< ?3 y+---+<iky N 


(5-7) 
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or in matrix form: 


S U, y) = { | x y ... 



<h 


(5-8) 


where polynomial terms are picked based on Pascal’s triangle in Table 5. 1. In the present 
capability polynomial order can range from 2 to 5. 

For least squares fitting, stresses in each CST are taken at the centroid of the element. 
Thus, for each CST element ‘i\ ‘x,’ and ‘yj’ refer to the element’s centroid position. Writ- 
ing polynomial equations for the stress a in ‘k’ elements leads to ‘k’ equations of the form 
[A s ]{q} = {b s } where 


1 x i y\ • 

■yi 

1 x 2 y 2 . 

A 

1 *3 >3 • 

-A 

_1 x k y k . 

-A 


(5-9) 


and 
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Table 5.1 - Choice of stress smoothing polynomial 
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{M = 


°K 


( 5 - 10 ) 


To solve for {q}, the normal equations approach is used (Ref. 18) to yield 


[A s ) T [A s ] { q } = [A s ) T {b s } 


(5-11) 


or 


[A 

new ] {<?} = {b 

new J 


(5-12) 


which can be directly solved using Ref. 16. 


5.4.2 Stress smoothing sensitivities 


Differentiating the previous equations for smoothed stresses with respect to a shape 
design variable |3 leads to 


dS (x, y) 


= { l,x,y, ...} 




dx dy ^ 

{0 ’5p’5fr" } {<?} 


( 5 - 13 ) 


where the vector { ( dq ) / (3|5) } is obtained from 



32 


, , 3 { it } ^[A new ] 

K 1 "3jr = — 3P — “ — 5p — {<,} 


Now 
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1 new J 
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(5-14) 


(5-15) 


(5-16) 


(5-17) 


and 
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9{/? 5 } 

^r 



W 


i ... 



( 5 - 18 ) 


For shape sensitivities, the expressions above take care of the motion of the (xj,yj) points 
used for least squares fitting as well as the motion of the point where the stress is calcu- 
lated. For sizing sensitivities, all points used for least squares fitting and stress output cal- 
culations are fixed and their derivatives are zero. Therefore, all derivatives of [A s ] with 
respect to any size variable tc are zero, resulting in 


[A new ] 


d{«} 

Bk 


3 k 


where 


( 5 - 19 ) 


_ r A i T ^A 

5tc 3 k 


( 5 - 20 ) 


and 


~3ir 



3 k 


( 5 - 21 ) 



CHAPTER 6 


WING MODEL TEST CASES 


6.1 Introduction 

Three wing models are described here for later use as test cases. For each, a brief phys- 
ical description is given along with all load cases to be examined. 

6.2 Gallagher wing 

The Gallagher model 1 wing (Ref. 19) is an unswept, untapered cantilever wing as 
shown in Figure 6.1. The aspect ratio is 4 and the depth-to-chord ratio is 0.075. All inter- 
nal members, being formed channels, are modeled as shear webs using membrane ele- 
ments. The channel flanges are modeled as rod elements whose cross-sectional area 
matches that of the flange area. Additionally, the skins arc modeled with membrane ele- 
ments only. The material used is 6061-T6 aluminum and its properties are: 

E= 10.0 x 10 ^ psi v = 0.3 p = 0.000259 lbm/in^ 

The load case analyzed is a 100 lbf. point load at each rib / spar intersection, first applied 
one at a time to derive the wings influence coefficients and then applied simultaneously 
(representing a continuous load over the wing) to examine its deformed shape. All wing 
skin thicknesses are 0.063” , web thicknesses are 0.040” and cap areas are modelled as 
being 0.02 square inches. Numerical tests include evaluation of the difference between 
modeling the spar / rib webs as plane stress elements carrying normal and shear stresses 
and spar / rib webs modeled by shear webs only. The effect of mesh refinement is exam- 
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Figure 6. 1 - Gallagher model 1 



- physical description 
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ined. The results from finite element models based on CST membranes are compared to 
those with LST membranes. 

6.3 Denke wing 

The Denke wing (Ref. 20) is a 45 degree swept back wing with an aspect ratio of 10. a 
depth-to-chord ratio of 0.35 and can be seen in Figure 6.2. Only four internal ribs are 
present along with the front and rear spars. Two load cases are considered. Load case 1 
involves a 1 lbf. point load applied vertically at the tip trailing edge, while load case 2 
involves a 1 lbf. point load applied vertically at the leading edge at 60% span. The mate- 
rial properties used are: 

E= 10.0 x 10 6 psi v = 0.3 p = 0.000259 lbm/in 3 

All wing skin thicknesses are 0.032”, web thicknesses are 0.051” and stringer areas are 
0.37 1 square inches for the leading and trailing edge elements, and 0.06 1 square inches for 
all remaining stringers. 

Again, effects of using plane stress and pure shear elements for spar and rib webs are 
studied as well as comparisons between the performance of CST’s and LST’s. This wing is 
an example of a thick, high aspect ratio wing typical in transonic transport airplane con- 
struction. Displacements and experimentally measured stresses in spar caps are used to 
assess accuracy of the present capability. 

6.4 TUmer/Martin/Weikel wing 

The Turner wing (Ref. 21), originally studied by Eggwertz and Noton, can be seen 
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physical description 
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in Figure 6.3. It has a 30 degree sweep, aspect ratio of 5 and a depth-to-chord ratio of 0.2 1. 
Five spars and three ribs are assumed to be perfectly attached to the upper and lower wing- 
skins and to each other. Cover skins are modeled as plane stress elements and a compari- 
son is made between modeling the spar and rib webs as plane stress or pure shear 
elements. The material used is aluminum with the following properties: 


E= 10.0 x 10 6 psi v = 0.3 p = 0.000259 Ibm/in 3 


All wing skin thicknesses are 0. 1 1 8”, web thicknesses are 0.059” and cap areas are 
0.0619 square inches. 

Measured displacements and skin stresses in the root area are used for evaluation. It 
should be mentioned (Ref. 21) that while measured skin stresses Q along the span (in 
the direction of the spars) are quite accurate, there is a reason to believe that normal 
stresses perpendicular to the spars ® xx and skin shear stresses 0 are inaccurate. Since 
they are very small compared to 0 , there would be no significant effect on failure esti- 

mation for the wing. 

For the Turner wing, in addition to the experimental data, finite element results, in par- 
ticular wing skin stresses and model natural frequencies, from a commercially available 
code (ELFINI, Ref. 22), were generated and used to compare to results from the present 
capability. 
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Figure 6.3 - Turner wing - physical description 




CHAPTER 7 


NUMERICAL RESULTS 


7.1 Introduction 


Three wing models are used to assess the present capability. First, accuracy of the 
finite element results needs to be evaluated, since the present capability is based on very 
basic, low order elements (in an effort to gain computational speed and obtain analytic 
sensitivities). This is done by comparing results obtained by the present capability to 
results by commercially available codes and to experimental results wherever possible. 

7.2 Gallagher wing 

Figure 7. 1 shows both the original wing skin mesh (based on existing ribs and spars) 
and a refined wing skin mesh employing four dummy ribs between each primary rib. 
When CSTs are used for cover skins and rib / spar webs, the effect of increasing the num- 
ber of spanwise divisions on the tip displacement using shear webs versus regular CST’s 
in the vertical webs is shown in Figure 7.2. Natural frequency convergence under mesh 
refinement is seen in Figure 7.3. As the number of divisions increase, the finite element 
displacement solution approaches that found experimentally (Ref. 1 9). It is interesting to 
note that the effect of modeling the spar and rib webs with shear webs becomes more 
important as the mesh is refined. The refined finite element model with shear webs is about 
5 % stiffer than the experimental model. 

A comparison of a refined CST model prediction (adrib = 5) and the LST model is 
shown in Figure 7.4. The LST model.uses a mesh based on the existing ribs and spars as in 
the coarser CST model. Mid-chord deflections along the entire span for both models are 
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Figure 7.1 - Gallagher model 1 wing skin meshes 
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Figure 7.2 - Mid-spar tip deflection convergence with span wise 
mesh refinement - Gallagher wing 
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Number of dummy ribs per section Number of dummy ribs per section 



Number of dummy ribc per section 


Figure 7.3 - Natural frequency convergence with spanwise 
mesh refinement - Gallagher wing 
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y position (in.) 


Figure 7.4 - Mid-spar deflection under a uniform load - Gallagher wing 
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compared to experimental data. Both models are in close agreement with only a 5.07% 
and a 2.05% wingtip deflection deviation from experiment, respectively (Table 7.1). 

Gallagher’s experimental influence coefficient matrix (Ref. 19) is reproduced in Table 

7.2 along with the resulting approximate influence matrix for the refined CST model 
(adrib=5) and the LST mode|. It can be seen that displacement results are good for the 
LST model while the CST model is slightly stiffen 

No experimental data is available for stresses on the Gallagher wing. As expected 
stresses in skin CST’s fluctuate and change discontinuously from element to element. Per- 
formance of the stress smoothing technique (Chapter 5) was evaluated by using polynomi- 
als of order two through five along with the adrib=4 mesh. The resulting polynomial fits 
are presented in Table 7.3, and plots along cuts A and B (Figure 7.1) are shown for each 
stress in Figures 7.5 through 7. 1 0. It is found that a polynomial of order N=4 captures CST 
stress variations well over the wing in this case. 

7.3 Denke wing 

Figure 7. 1 1 shows both the original wing skin mesh (based on existing spars and ribs) 
and a refined wing skin mesh employing two dummy ribs between each pair of primary 
ribs. Deflection results for the Denke wing (in the case of CST elements) with an increas- 
ing number of spanwise divisions are compared in Figures 7.12 and 7.13. Results of using 
shear webs and CST membranes (including normal stresses) for spar and rib webs are 
shown for both load cases. Natural frequency convergence results are shown in Figure 
7.14. 

Excellent correlation between experiment and finite element modeling using CST’s is 
shown in Figures 7.15 and 7.16 for load cases 1 and 2, respectively. The CST model used 
for these and all subsequent results has adrib=2. Comparison of results from the LST 
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Table 7. 1 - Displacements of the Gallagher model 1 wing 


Mid-chord deflection (in.) 
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CST model 

% error 

LST model 
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3.0 

1.319 

1.252 
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2.05 

8.0 
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9.65 
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3.24 

13.0 
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14.17 

0.233 

9.51 
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Tabic 7.2 - Gallagher model 1 influence coefficients 
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0.0262 





0.1*63 

0.06*3 

0.0764 

0.0*4 

0.0913 

0.090 

0.0169 

0X3213 

00247 

0.0274 

0.0293 






0.0631 

o.osa 

0.0473 

0.0411 

0X335 

0X3215 

0.0115 

0.0153 

0.0123 

0.009* 







0.0521 

0.04*6 

0.044* 

0.0411 

0.01* 

0X3172 

0X3156 

0.0139 

0.0116 








0.0497 

0.04*6 

0.0475 

0X1147 

0.0155 

0.015* 

0.0155 

0.0146 









0.052* 

0.054* 

0X1126 

0X3139 

0.0156 

0.0172 

0.01* 










0X3631 

0X101* 

0X3123 

0X3153 

0.01*5 

0.0215 

10 










0j0125 

0X30*6 

0X305* 

0X303* 

0.002 

II 
• <* 











0X10*5 

0X1065 

0X3051 

0X303* 


0.0074 0.0065 0.0051 

0.0015 0.0016 

0.0125 


LST model 
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Table 7.3 - Gallagher model 1 stress smoothing polynomials 


Sigma X: 

N=2 S = 4144 + 163. lx - 434.5y - 1 1 ,72x 2 + 0.85xy + 9.85y 2 

N=3 S = 4600 + 523. 8x - 870y - 38. 9x 2 - 17.7xy +49.1y 2 + 0.145x 3 + 1.51x 2 y- 

0. 137xy 2 - 0.85y 3 

N=4 S = 6279 - 1 16. 6x - 1564y + 130.8x 2 - 20.62xy + 153. 8y 2 - 15.35x 3 + 
0.382x 2 y + 0. 195xy 2 - 6.255y 3 + 0.447x 4 + 0.139x 3 y - 0.067x 2 y 2 + 
0.0149xy 3 + 0.088y 4 

N=5 S = 8613 - 968.9x - 2890y + 229.6x 2 + 168.9xy + 406.5y 2 -5.66x 3 - 

34. 1 5x 2 y - 6.05xy 2 - 27.6y 3 * 1 .52x 4 + 2.56x 3 y + 0.6x 2 y 2 + 0. 1 1 6xy 3 + 
0.87y 4 + 0.075x 5 - 0.06x 4 y - 0.025x 3 y 2 - 0.0024x 2 y 3 - 0.00 lxy 4 - 
0.0 ly 5 

Sigma Y: 

N=2 S = 25500 + 54.36x - 1478y - 2.72x 2 - 0.43xy + 21.2y 2 

N=3 S = 24890 + 396.3x - 1467y - 44.6x 2 - 10.4xy + 22.5y 2 + 1.45x 3 + 0.61x 2 y + 

O^xy 2 - 0.034y 3 

N=4 S = 24460 + 323. 4x - 1245y + 18.06x 2 - 51.6xy + 2.16y 2 - 6.38x 3 + 

OJSx^ + 1.53xy 2 + 0.7 89y 3 + 0.27 lx 4 -0.02x 3 y - 0.076x 2 y 2 - 
O.OOSxy 3 - 0.0 13y 4 

N=5 S = 24160 + 998. 5x - 1334y - 389.7x 2 + 37xy + 4.16y 2 + 81.75x 3 - 
20.4x 2 y + 1.92xy 2 + 0.387y 3 - 7.39x 4 + 2.13x 3 y + 0.065x 2 y 2 - 
0.072xy 3 + 0.0 12y 4 + 0.23x 5 - 0.066x 4 y - 0.0062x 3 y 2 - 0.0x 2 y 3 + 

0.001 lxy 4 - 0.00044y 5 

Tau XY: 

N=2 S = -1 149 + 144.2x + 86.33y - 0.0584x 2 - 10.45xy - O.By 2 

N=3 S = -2435 + 36.5x + 446y + 54.9x 2 - 56.4xy - By 2 - 2.4x 3 + 0x 2 y + 
l^xy 2 - 0.03y 3 

N=4 S = -3819 + 78.6x + 1085y + 94.8x 2 - 134.2xy - 69.7y 2 - 6.3x 3 - O.^x^ + 
8-Bxy 2 + 1 .5y 3 + 0. 1 28x 4 + 0.0026x 3 y + 0.0043x 2 y 2 - 0. 148xy 3 - 
0.006y 4 

N=5 S = -6531 + 1865x + 2029y - 473.4x 2 - 309.5xy - 207y 2 + 84. 1 x 3 + 9.8x 2 y + 
26. lxy 2 + 9.27^ - 6.66x 4 - 0.0013x 3 y - 0.577xV - 0.9xy 3 - 0. 18y 4 + 

0.1 9x 5 + 0.028x 4 y + 0.028x 3 y 2 - 0.0012x 2 y 3 + 0.013xy 4 + 0.00 ly 5 
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Figure 7.5 - O w stress smoothing along line A - Gallagher wing 
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Figure 7.6 - <7 yy stress smoothing along line A - Gallagher wing 
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Figure 7.7 - a xy stress 
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Figure 7.9 - Oyy stress smoothing along line B - Gallagher wing 
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Figure 7. 10 - C xy stress smoothing along line B - Gallagher wing 
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Figure 7.12 - Trailing edge tip deflection convergence with spanwise 
mesh refinement - Denke wing (load case 1) 
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Figure 7.13 - Trailing edge tip deflection convergence with spanwise 
mesh refinement - Denke wing (load case 2) 
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Figure 7.14- Natural frequency convergence with spanwise 
mesh refinement - Denke wing 
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Figure 7.15 - CST model leading and trailing edge deflection - Denke wing (load case 1) 
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Figure 7. 16 - CST model leading and trailing edge deflection - Denke wing (load case 2) 
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model (basic mesh) with experimental results are shown in Figures 7.17 and 7. IS. The 
correlation is good. It is slightly stiffer than the best case CST model results along the 
leading edge in load case 2 (Figure 7. 16). Numerical results are listed in Table 7.4 and are 
seen to be reasonably close to experiment. 

Stress behavior is first analyzed in the spar caps along the wing root. Figure 7.19 
shows the comparison for each model with respect to published values (Ret. 2(1). As can 
be seen, finite element stress magnitudes are lower than experiment towards the root trail- 
ing edge, and a rather large discrepancy exists in the CST model at 80% ot the chord. One 
reason for this is that the Ref. 20 results were taken at the wing root while the finite ele- 
ment results were taken from the mid-point of the root rod element (the axial stress is 
assumed constant throughout its length). In addition, the well known root trailing edge 
stress singularity in swept back wings appears, and accuracy of all calculated results dete- 
riorates in that region. 

A comparison of cap stresses along the leading and trailing edge spar caps for each 
model (CST and LST) for load case 2 is seen in Figures 7.20 and 7.21. Good correlation 
with experiment exists. The stress values plotted for each spar cap element were taken at 
its geometric midpoint as mentioned previously. 

No experimental wing skin stress data is available for the Denke wing. Skin stress 
curve fits were again attempted for each stress along both a chordwise and a spanwise cut 
(Fig. 7.1 1) in the CST model. Numerical details of the various curve fitting polynomials 
are given in Table 7.5. Figures 7.22 through 7.27 show the resulting plots. As with the 
Gallagher wing, a polynomial of order N=4 gives the best representation for the fluctuat- 
ing CST element stresses with Gyy showing the best behavior. 
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Figure 7.17 - LST model leading and trailing edge deflection - Denkc wing (load case 1) 
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Figure 7.18 - LST model leading and trailing edge deflection - Dcnke wing (load case 2) 
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Table 7.4 - Displacements of the Denke wing 


Load case 1 Deflection (xlO-3 in.) 


Node 

experiment 

CST model 

% error 

LST model 

% error 

1 



0.00 

0.007 

0.00 

2 


0.035 


0.035 

6.06 

3 

0.040 

0.040 

0.00 

0.039 

2.50 

4 

0.096 

0.102 

6.25 

0.099 

3.13 

5 

0.100 

0.097 

3.00 

0.094 

6.00 

6 

0.185 

0.189 

2.16 

0.183 

1.08 

7 

0.170 

0.174 

2.35 

0.168 

1.18 

8 

0.290 

0.291 

0.34 

0.281 

3.10 

9 

0.260 

0.263 

1.15 

0.255 

1.92 

10 

0.400 

0.403 

0.75 

0.389 

2.75 

Load case 2 
Node 

experiment 

Deflection (xlO-' 
CST model % error 

) in.) 

LST model 

% error 

1 

0.014 


7.14 

0.012 

14.29 

2 

0.011 

1 v' 

18.18 

0.012 

9.09 

3 

0.038 

I 


0.036 

5.26 

4 

0.032 



0.033 

3.13 

5 

0.070 

0.067 

4.29 

0.066 

5.71 

6 

0.056 

0.056 

0.00 

0.054 

3.57 

7 

0.089 

0.087 

2.25 

0.084 

5.62 

8 

0.076 

0.076 

0.00 

0.074 

2.63 

9 

0.108 

0.107 

0.93 

0.104 

3.70 

10 

0.096 

0.097 

1.04 

0.094 

2.08 





axial stress (psi) 
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Figure 7.19- Root chord cap stresses - Denke wing (load case 2) 
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Figure 7.20 - CST model spar cap stresses - Denke wing (load case 2) 
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Figure 7.2 1 - LST model spar cap stresses - Denke wing (load case 2) 
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Table 7.5 - Denke stress smoothing polynomials 


Sigma X: 

N=2 S = -0.5732 + O.I33x - 0.2246y - 0.0005606x 2 - 0.00098 17xy + 0.004 109y 2 

N=3 S = 0.9 13 1 - 0.05502x - 0.2869y + 0.007494x 2 - 0.00063xy + 0.005896y 2 - 

0.000 1728x 3 + 0.000346x 2 y - 0.0004927xy 2 + 0.0002019^ 

N=4 S = 1.616 - 0.259 lx - 0.2 163y + 0.02552x 2 - 0.00929xy + 0.002 164^ - 
0.0009427x 3 + 0.0016 16x 2 y - 0.00 18xy 2 + 0.000973}^ + 0.00001 1 7x 4 - 
0.00002935x 3 y + 0.0000317x 2 y 2 - 0. 000009087 xy 3 - 0.0000059ily 4 
N=5 S = 2. 1 27 - 0.98x + 0.5835y + 0.25 1 x 2 - 0.467xy + 0.228y 2 - 0.028x 3 + 

0.078 16x 2 y - 0.07174xy 2 + 0.02 13y 3 + 0.0014x 4 - 0.005 123x 3 y + 
0.00692x 2 y 2 - 0.004 126xy 3 + 0.000937y 4 - 0.0000258x 5 + 0.0001 15x 4 y - 
0.0002024x 3 y 2 + 0.000 175x 2 y 3 - 0.000074xy 4 + 0.00001 185y 5 

Sigma Y: 

N=2 S = 4.56 - 0.0895x - 0.08 107y - 0.0004309x 2 + 0.001 836xy + 0.0001947)^ 

N=3 S = 2.713 + 0.1 17 lx + 0.0265y - 0.008532x 2 - 0.003825xy + 0.002099y 2 + 

0.0003573x 3 - 0.001 1 19x 2 y + 0.00 Hxy 2 - 0.000823 ly 3 
N=4 S = 1.683 + 0.1255x + 0.3875y + 0.00897x 2 - 0.06642xy + 0.02704y 2 - 

0.00097 x 3 + 0.002349x 2 y + 0.0007279xy 2 - 0.001396)^ + 0.00004581x 4 - 
0.000 194x 3 y + 0.0003256x 2 y 2 - 0.0002886xy 3 + 0.0001078y 4 
N=5 S = 2.161 - 0.098x + 0.539y - 0.0287x 2 + 0.077xy - 0.0864JT 2 + 0.01676x 3 - 
0.06263x 2 y + 0.07636xy 2 - 0.02968y 3 - 0.0017 12x 4 + 0.00766x 3 y - 
0.01253x 2 y 2 + 0.008884xy 3 - 0.00230 ly 4 + 0.0000539x 5 - 0.0002943x 4 y + 
0.0006389x 3 y 2 - 0.0006924xV+ 0.0003766xy 4 - 0.00008273^ 

Tau XY: 

N=2 S = -1.191 + 0.00433x + 0.05274y - 0.0002078x 2 + 0.000339xy - 0.000801y 2 

N=3 S = -1.368 + 0.1041x - 0.07157y - 0.004549x 2 + 0.002929xy + 0.00359^ - 

0.0000683x 3 + 0.0004915x 2 y - 0.0007492xy 2 + 0.0002806)^ 

N=4 S = -1.556 + 0.1912x - 0.1758y - 0.002212X 2 - 0.02199xy + 0.03325^ - 

0.0006887x 3 + 0.003224x 2 y - 0.003743xy 2 + 0.0008293y 3 + 0.0000 1036x 4 - 
0.00002879x 3 y - 0.000003943x 2 y 2 + 0.00005072xy 3 - 0.00002422y 4 
N=5 S = - 1 .63 + 0.086x + 0.03225y + 0.0284x 2 - 0.08423xy + 0.05 1 24^ - 
0.0043x 3 + 0.0126x 2 y - 0.01077xy 2 + 0.002693^ + 0.0002488x 4 - 
0.00095x 3 y + 0.00137x 2 y 2 - 0.0009256xy 3 + 0.00025y 4 - 0.0000067x 5 + 
0.0000356x 4 y - 0.000078x 3 y 2 + 0.00009xV - 0.00005xy 4 + 0.00001 ly 5 




Figure 7.22 - O u stress smoothing along line A - Denke wing (load case 2) 
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Figure 7.23 - <Jyy stress smoothing along line A - Denke wing (load case 2) 
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Figure 7.25 - O xx stress smoothing along line B - Denke wing (load case 2) 
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Figure 7.26 - Oyy stress smoothing along line B - Denke wing (load case 2) 
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7.4 Ttirner/Martin/Weikel wing 


Figure 7.28 shows both the basic wing skin mesh and the refined wing skin mesh for 
the Turner wing using CST elements. In this case the mesh includes both spanwise and 
chordwise refinements, introducing the presence of floating nodes. Chordwise mesh 
refinement consists of one dummy spar per spar interval, and is employed to allow for 
more CST elements across the chord. One dummy rib per rib interval is then added within 
the root region. Figure 7.29 shows the LST element model used for comparison. 

The effect on wing deflection of modeling the Turner wing with a refined mesh and 
pure CST shear webs as compared to the LST model is seen in Figures 7.30 and 7.31. 
Mesh refinement has only a small effect on the spanwise vertical deflection. Tables 7.6 and 
7.7, though, show that refining the mesh in this case leads to greater in-plane deflections 
(x- and y- axes) for the CST model. 

With respect to stresses, nodal stress averaging following the results of Turner (Ref. 
21) was performed for each CST wing model. Tables 7.8 and 7.9 contain the nodal aver- 
ages for each of the three stresses within the wing root area as compared to published 
results. Close agreement is found for <J XX and G yy . In the case of the shear stress O xy the 
correlation is not as good. 

Natural frequency results for both finite element models as compared with those avail- 
able from a commercial finite element package (ELF1NI, Ref. 22) can be seen in Table 
7.10. Excellent agreement using the original mesh is evident Natural frequency results 
using the refined mesh decrease in accuracy as the frequency increases due to localized 
vibration of lumped masses at floating nodes. An attempt to solve this problem involved 
studying the choice of dummy element thicknesses (1% of a real element’s thickness was 
the choice in all studies up to this point). The effect of varying the dummy element thick- 
ness from 1% to 10% on displacements and natural frequencies is shown in Table 7.11. 



Final CST wingskin raesh (Adspar=l) 


Figure 7.28 - Turner wing skin CST meshes 








Figure 7.30 - Leading edge deflection - Turner wing 
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Tabic 7.6 - Displacements of the Turner wing (original mesh) 


Turner (x 10-6) 


CST model (xlO-6) 


0.008 

-0.443 

-0.850 

-1.225 

-1.585 

-0.030 

-0.660 

-1.014 

-1.327 

-1.682 

0 

-0.480 

-0.858 

-1.184 

-1.511 

0 

-0.513 

-0.904 

-1.226 

0 

-0.483 

-0.836 

0 

-0.383 

0 


-4.491 

-4.333 

-4.251 

-4.142 

-4.060 

- 2.666 

-2.840 

-3.043 

-3.069 

-2.991 

0 

-1.252 

-1.801 

-2.029 

-2.098 

0 

-0.940 

-1.341 

-1.548 

0 

-0.679 

-1.016 

0 

-0.506 

0 


-15.910 

-16.690 

-16.939 

-16.069 

-13.669 

-5.695 

-7.797 

-8.947 

-8.933 

-7.463 

0 

-2.004 

-3.421 

-3.923 

-3.232 

0 

-1.091 

-1.689 

-1.333 

0 

-0.384 

-0.178 

0 

0.260 

0 


.009 
1.422 
1.820 
-1.180 
-1.532 
0.009 
-0.616 
-0.983 
-1.299 
-1.662 
0 

-0.393 

-0.798 

-1.141 

-1.478 

0 

-0.440 

-0.850 

-1.191 

0 

-0.429 

-0.788 

0 

-0.338 

0 


-4.874 

-4.747 

-4.650 

-4.530 

-4.473 

-2.878 

-3.104 

-3.329 

-3.368 

-3.311 

0 

-1.315 

-1.960 

-2.231 

-2.328 

0 

-0.981 

-1.455 

- 1.688 

0 

-0.708 

-1.076 

0 

-0.504 

0 
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Tabic 7.7 - Displacements of the Turner wing (refined mesh) 


Node 

Turner (xl0-6) 

CST model (xl0-6) 

u 

V 

w 

u 

V 

w 

31 

0.008 

-4.491 

-15.910 

0.103 

-5.074 

-18.053 

32 

-0.443 

-4.333 

-16.690 

-0.425 

-4.857 

-19.000 

33 

-0.850 

-4.251 

-16.939 

-0.828 

- 4.744 

-19.312 

34 

-1.225 

-4.142 

-16.069 

-1.198 

-4.616 

-18.458 

35 

-1.585 

-4.060 

-13.669 

-1.562 

-4.562 

-16.005 

36 

-0.030 

-2.666 

-5.695 

0.048 

-3.066 

-6.520 

37 

-0.660 

-2.840 

-7.797 

-0.629 

-3.227 

-9.634 

38 

-1.014 

-3.043 

-8.947 

-0.993 

-3.410 

-10.190 

39 

-1.327 

-3.069 

-8.933 

-1.304 

-3.429 

-10.259 

40 

-1.682 

-2.991 

-7.463 

-1.658 

-3.363 

-8.665 

41 

0 

0 

0 

0 

0 

0 

42 

-0.480 

-1.252 

-2.004 

-0.444 

-1.426 

-2.133 

43 

-0.858 

-1.801 

-3.421 

-0.821 

-2.054 

-3.772 

44 

-1.184 

-2.029 

-3.923 

-1.144 

-2.300 

-4.400 

45 

-1.511 

-2.098 

-3.232 

-1.476 

-2.379 

-3.588 

46 

0 

0 

0 

0 

0 

0 

47 

-0.513 

-0.940 

-1.091 

-0.485 

-1.039 

-1.196 

48 

-0.904 

-1.341 

-1.689 

-0.874 

-1.502 

-1.877 

49 

-1.226 

-1.548 

-1.333 

-1.208 

-1.713 

-1.398 

50 

0 

0 

0 

0 

0 

0 

51 

-0.483 

-0.679 

-0.384 

-0.467 

-0.733 

-0.406 

52 

-0.836 

-1.016 

-0.178 

-0.822 

-1.089 

-0.089 

53 

0 

0 

0 

0 

0 

0 

54 

-0.383 

-0.506 

0.260 

-0.377 

-0.503 

0.349 

55 

i 

0 

0 

0 

0 

0 

0 
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Tabic 7.8 - Summary of Turner computed nodal stress averages (original mesh) 


Turner stress averages (psi) 


CST stress averages (psi) 


Sigma X Sigma Y Sigma XY 


Sigma X Sigma Y Sigma XY 


Node 


36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 


0.36 

0.61 

0.60 

0.34 

0.20 

0.98 

0.91 

0.59 

0.29 

0.15 

0.52 

0.29 

0.15 

0.08 

-0.07 

-0.08 

-0.05 

-0.26 

- 0.22 


8.93 

6.68 

5.41 

4.60 
4.28 

9.55 
7.43 

5.56 

4.42 
3.80 
7.13 

5.60 
4.27 
3.55 
5.18 
4.10 

3.43 
3.71 
3.25 


0.45 

-0.38 

-0.17 

0.04 

0.12 

-0.87 

-0.41 

0.00 

0.25 

0.38 

-0.15 

0.17 

0.41 

0.52 

0.24 

0.44 

0.53 

0.37 

0.41 


0.49 

0.61 

0.51 

0.30 

0.13 

1.05 

0.90 

0.51 

0.21 

0.09 

0.66 

0.29 

0.10 

0.00 

0.06 

-0.09 

-0.09 

-0.16 

-0.24 


8.85 
7.57 
5.73 

4.85 

4.42 

9.11 

7.42 
5.94 
4.71 

4.11 
6.92 
5.64 
4.59 
3.91 
5.03 
4.11 
3.67 
3.61 
3.02 


-0.28 

-0.51 

-0.19 

- 0.01 

0.12 

- 0.88 

-0.58 

-0.07 

0.23 

0.34 

-0.39 

-0.03 

0.33 

0.47 

0.04 

0.28 

0.43 

0.25 

0.19 
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Tabic 7.9 - Summary of Turner computed nodal stress averages (refined mesh) 
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Tabic 7.10 - Natural frequencies of the Turner wing 


Natural Frequency (Hz) 


Mode 

Original 
CST model 

Refined 
CST model 

ELFINI 

1 

120 

119 

116 

2 

337 

327 

318 

3 

419 

411 

418 

4 

602 

540 

577 

5 

1107 

687 

1086 
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Tabic 7. 1 1 - Dummy thickness effect on Turner displacements and 

natural frequencies 


Natural Frequencies (Hz) 


Mode 

Original 
CST model 

Refined CST models 
(1%) (5%) (10%) 

ELFINl 

1 

120 

119 

119 

119 

116 

2 

337 

327 

327 

326 

318 

3 

419 

411 

426 

429 

418 

4 

602 

540 

593 

599 

577 

5 

1107 

687 

1074 

1130 

1076 

Load point 
displacement 
(10* 3 in.) 

0.216 

0.219 

0.218 

0.217 

N/A 









Increased thicknesses result in increased natural frequency accuracy with an insignificant 
decrease in tip displacement. The study’s effect on element stresses can be seen in Figures 
7.32 and 7.33. Figure 7.32 details the change in leading and trailing edge cap stresses for 
an increase in dummy membrane thickness while Figure 7.33 shows the change in CST 
element stresses and G yy along line B 2 . Again, increased dummy element thicknesses 
have a minimal effect on both cap and membrane stresses. Subsequently, a good rule of 
thumb is to use dummy elements with a thickness of between 5% and 10% of what the 
actual structure requires only if accurate natural frequency / mode shape information is 
desired. 

Stress smoothing was again employed for each CST wing model, with finite element 
stress results available from ELFIN1. For each stress, a polynomial was found at each cut 
A and Bj/E^ (Figures 7.28 and 7.29) for both the basic CST model and the LST model 
while being compared to ELFINI results and the CST model’s smoothed stresses. Numer- 
ical details of the various curve fitting polynomials for each wing mesh are in Tables 7. 1 2 
and 7. 13, but with a polynomial order of N=4 having been previously established, this 
degree will be used for all subsequent stress comparisons. Plots are shown in Figures 7.34 
through 7.39. Note the linear, piecewise continuous nature of the LST model’s stresses 
along each cut. Additionally, the ELFINI stress results can be seen in Figures 7.40 through 
7.42 for each stress. 

Good agreement with ELFINI for all models can easily be seen for Gyy along with 
excellent curve fits at each cut. Along cut A, reasonable accuracy in both the element 
stresses and the N=4 stress polynomial is obtained for G xx while poor accuracy between 
ELFINI and stress smoothing results exist for G xy . It is also worth noting the decrease in 
accuracy as one nears the trailing edge root location (100% chord). Along spanwise cuts 
Bj/B 2, the N=4 curve fit and element stresses are better for G xy but this time G M ELFINI 
results and stress results are quite different. In general, good agreement exists between the 
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Figure 7.32 - Dummy thickness effect on Turner cap stresses 
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Figure 7.33 - Dummy thickness effect on Turner membrane stresses along line B2 
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Tabic 7. 12 - Turner stress smoothing polynomials (original CST mesh) 


Sigma X: 

N=2 S = 0.04528 + 0.04588x - 0.003832y - 0.00 19x 2 - 0.003465xy + 0.0009846)^ 

N=3 S = 0.6845 - 0.03744x - 0. 1 5 1 1 y + 0.00983x 2 + 0.004395xy + 0.008456y 2 - 

0.001 147x 3 + 0.0006225x 2 y - 0.0004522xy 2 - 0.00008472y 3 
N=4 S = 0.908 - 0.4852x - 0.07524y + 0. 1 462x 2 + 0.0 1 168xy - 0.00089 1 y 2 - 

0.01702x 3 - 0.001019x 2 y - 0.0001546xy 2 + 0.0002379^ + 0.0006553x 4 + 
0.0x 3 y + 0.00005336x 2 y 2 - 0.0000 195xy 3 - 0.000002632y 4 
N=5 S = 1 .534 - 1 ,709x - 0.08587y + 0.9 1 85x 2 - 0. 1 164xy + 0.02993y 2 - 0. 1907x 3 + 
0.01492x 2 y + 0.007555xy 2 - 0.002937y 3 + 0.017 19x 4 - 0.000755x 3 y - 
0.000653x 2 y 2 - 0.000 ISOSxy 3 + 0.000 107y 4 - 0.00057 16x 5 + 0.000033x 4 y - 
0.000004953x 3 y 2 + 0.0000 1668x 2 y 3 - 0.00000 1182xy 4 - 0.00000 1176y 5 

Sigma Y: 

N=2 S = -1.1 +0.01 66x + 0.41 74y + 0.01 546x 2 - 0.0238xy - 0.002058y 2 
N=3 S = 0.6766 - 0. 182 lx - 0.0245 ly + 0.0606 lx 2 - 0.007962xy + 0.02102^ - 
0.004979x 3 + 0.002886x 2 y - 0.001505xy 2 - 0.0002467y 3 
N=4 S = 0.9147 - 0.979x + 0.207 ly + 0.3286x 2 + 0.01566xy - 0.01273y 2 - 
0.044x 3 + 0.004483x 2 y - 0.002706xy 2 + 0.001263^ + 0.001983x 4 - 
0.0005238x 3 y + 0.0002468x 2 y 2 - 0.00003682xy 3 - 0.0000 182y 4 
N=5 S= 1.854 - 3.575x + 0.5373y + 1.824x 2 - 0.1196xy- 0.0328 ly 2 -0.388x 3 + 
0.023 12x 2 y + 0.00878xy 2 + 0.0005012y 3 + 0.03721x 4 - 0.004233x 3 y + 
0.000665x 2 y 2 - 0.00059 Uxy 3 + 0.00006662y 4 - 0.00132x 5 + 0.000253x 4 y - 
0.00007632x 3 y 2 + 0.0000 1977x 2 y 3 + 0.000004684xy 4 - 0.000001425)^ 

Tau XY: 

N=2 S = -1.592 + 0.07 176x + 0.0341y + 0.006499x 2 - 0.00861 lxy + 0.0017 13y 2 

N=3 S = -0.9605 + 0.0574x - 0. 1024y - 0.01467x 2 + 0.009038xy + 0.006647y 2 + 

0.001343x 3 - 0.00002344x 2 y - 0.0005 ^xy 2 - 0.00002641^ 

N=4 S = -1.064 - 0.04879x + 0.03493y - 0.03781x 2 + 0.03417xy - 0.01263)^ + 

0.0098 17x 3 - 0.005 126x 2 y - 0.0002034xy 2 + 0.0007 198y 3 - 0.0004908x 4 + 
0.000 1609x 3 y + 0.00007006x 2 y 2 - 0.00002309xy 3 - 0.00000845y 4 
N=5 S = -2.28 + 0.944x + 0.37 1 y - 0.423x 2 - 0.06 1 24xy - 0.0476^ + 0.08784x 3 - 
0.001 14x 2 y + O.OOSSSSxy 2 + 0.001904y 3 - 0.007789x 4 + 0.0003322x 3 y - 
0.000295x 2 y 2 - 0.0003 157xy 3 - 0.0000147y 4 + 0.00025x 5 - 0.00000504x 4 y - 
0.000003095x 3 y 2 + 0.000008086x 2 y 3 + 0.000002909xy 4 - 0.000000 1938y 5 



Table 7. 13 - Turner stress smoothing polynomials (refined CST mesh) 


Sigma X: 

N=2 S = -0.039 + 0.06141x + 0.002043y - 0.002725x 2 - 0.003556xy + 0.0008221y 2 
N=3 S = 0.5 1 1 1 + 0. 1 366x - 0. 1 68 1 y - 0.0227 1 x 2 + 0.003742xy + 0.00983/* + 
0.0006 10 lx 3 + 0.0005858x 2 y - 0.0004033xy 2 - 0.000 1182y 3 
N=4 S = 0.501 - 0.02049x - 0.07597y - 0.001607x 2 + 0.01968xy - 0.0035 ly 2 + 

0.000373x 3 - 0.001673x 2 y - 0.0005401xy 2 + 0.0004465/ - 0.00003432x 4 + 
0.00003975x 3 y + 0.00004602x 2 y 2 - 0.000007925xy 3 - 0.00000725 ly 4 
N=5 S = 0. 1 862 + 0. 1 0 1 1 x - 0.066 1 1 y + 0. 1 1 74x 2 - 0. 1 1 33xy + 0.02747y 2 - 
0.032 18x 3 + 0.008646x 2 y + 0.009438xy 2 - 0.003 ly 3 + 0.002985x 4 + 

0.000 123x 3 y - 0.0006824x 2 y 2 - 0.000238xy 3 + 0.000 1248y 4 - 0.0000922x 5 - 
0.0000079x 4 y + 0.0x 3 y 2 + 0.0000 145x 2 y 3 + 0.00000079xy 4 - 0.000001 59/ 

Sigma Y: 

N=2 S = -1.239 + 0.04716x + 0.4226y + 0.01422x 2 - 0.02482xy - 0.00193y 2 
N=3 S = 0.4089 + 0.1094x - 0.03 197y - 0.005225x 2 - 0.001642xy + 0.01977/ - 
0.00l089x 3 + 0.002496x 2 y - 0.001543xy 2 - 0.0002065y 3 
N=4 S = 0.3438 - 0.1477x + 0.1505y + 0.03614x 2 + 0.03247xy - 0.007506y 2 - 

0.005673x 3 - 0.00005743x 2 y - 0.002 1 26x/ + 0.0009 176/ + 0.0003089x 4 - 
0.0002046x 3 y + 0.000 1907x 2 y 2 - 0.00003542xy 3 - 0.0000 1229y 4 
N=5 S = 0.02015 - 0.626x + 0.4174y + 0.395x 2 - 0.04725xy - 0.0265y 2 - 0.088x 3 + 
0.0000423 lx 2 y + 0.008073xy 2 + 0.000368/ + 0.0086x 4 - 0.0006335x 3 y + 
0.000297x 2 y 2 - 0.00049xy 3 + 0.00005806y 4 - 0.000307x 5 + 0.0000548x 4 y - 
0.00002804x 3 y 2 + 0.0000075 1 9x 2 y 3 + 0.000005508xy 4 - 0.00000 1 306/ 

Tau XY: 

N=2 S = -1.782 + 0.1 22x + 0.0442y + 0.003449x 2 - 0.009177xy + 0.00144l/ 

N=3 S = - 1 .092 + 0. 1442x - 0. 108y - 0.02863x 2 + 0.009283xy + 0.007 144/* + 
0.002026x 3 - 0.00008864x 2 y - 0.0004913xy 2 - 0.00004331y 3 
N=4 S = -1.237 + 0.02735x - 0.02508y - 0.1311x 2 + 0.03314xy - 0.005875y 2 + 
0.01912x 3 - 0.003733x 2 y - 0.000557xy 2 + 0.0004753y 3 - 0.0007873x 4 + 
0.00007 442x 3 y + 0.00006826x 2 / - 0.00001512xy 3 - 0.0000058 12y 4 
N=5 S = -2.327 + 0.8856x + 0.4347y - 0.422x 2 - 0.042xy - 0.0623 ly 2 + 0.0937x 3 - 
0.009295x 2 y + 0.009 168xy 2 + 0.003004 y 3 - 0.008746x 4 + 0.00123x 3 y - 
0.0002 126x 2 y 2 - 0.0003394xy 3 - 0.0000504y 4 + 0.00029x 5 - 0.0000345x 4 y - 
0.0000 1082x 3 y 2 + 0.000008993x 2 y 3 + 0.000003068xy 4 + 0.0000002293/ 
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Original CST mesh 



Final CST mesh 
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Figure 7.34 - O u stress smoothing along line A - Turner wing 
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Original CST mesh 



Final CST mesh 
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Figure 7.35 - G yy stress smoothing along line A - Turner wing 
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Original CST mesh 



Final CST mesh 



% span 


Figure 7.37 - O u stress smoothing along lines Bj and B2 - Turner wing 
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Figure 7.38- Oy y stress smoothing along lines Bj and B2 - Turner wing 
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Original CST mesh 



Final CST mesh 
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Figure 7.39 - O xy stress smoothing along lines Bj and B 2 - Turner wing 
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Figure 7.41 - Oyy stress contour plot - ELFINI finite element model 
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Figure 7.42 - O xy stress contour plot - ELFINI finite element model 
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point stresses of the CST elements and the linear stresses of the LST elements. 

It should be remembered, however, that both O xx and G xy are significantly small as 
compared with G yy . Thus, failure predictions for the Turner wing by the current CST/LST 
modeling technique and ELFIN1 as well as test data will all be in good agreement. Also, 
there is a doubt as to the accuracy of measured G xx and G xy values, and large stress gradi- 
ents at the root trailing edge are certainly affecting accuracy of these small stresses. 



CHAPTER 8 


ANALYTIC SENSITIVITY RESULTS 


8.1 Introduction 


Analytic sensitivity calculations are checked by corresponding finite difference deriv- 
atives. In addition, computational efficiency issues of employing analytic sensitivities ver- 
sus finite difference sensitivities is evaluated. The wing models of choice for all future 
discussions are the Gallagher model 1 wing (adrib = 4) and the Denke wing (adrib=2) both 
having shear web CSTs. 

8.2 Analytic sensitivities vs. finite difference sensitivities 

With respect to finite difference methods, the expression 

dx Ax x->~x, 

if = _ = J i (8-1) 

ov Av v 2 - Vj 

describes the derivative of any behavior function *x’ with respect to a change in any vari- 
able ‘v.’ For large perturbations in *v,’ truncation error results in inaccurate derivatives 
due to it being a first order approximation, while theoretically as Av approaches zero, the 
approximation becomes exact. Realistically, this process introduces round-off errors due 
to computer finite length representation of numbers (Ref. 1). 

As an example of shape design variable sensitivity, the Gallagher model under a uni- 
form load and it’s perturbed version with respect to both xpR and y^ arc compared. The 
analytic sensitivities of the vertical displacement at the trailing edge tip, the first natural 
frequency, the leading edge root cap stress and the spanwise plane stress Oyy in the CST 
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leading edge wing skin root element are calculated. In using Hnite differences, perturba- 
tions of 0.001% to 0.1% of the characteristic dimension (chord length for xp-R and span 
length for y R ) are used. The results are found in Tables 8.1 and 8.2. The Denke wing 
model under a 100 lb. trailing edge tip load is tested in the exact same fashion as above. 
Results are shown in Tables 8.3 and 8.4. Since the program is written in double precision, 
round-off errors in the finite difference scheme for small perturbations do not show for the 
range analyzed. For larger perturbations, truncation error explains any discrepancies. The 
analytic sensitivities are seen to be in complete agreement. 

As an example of sizing design variable sensitivity, the Gallagher model is used and 
the same sensitivities are sought, this time with respect to the cross-sectional area of the 
leading edge spar cap element. Table 8.5 shows the results. Again, the same performance 
as detailed for the shape sensitivities is achieved. 

To further exhibit the accuracy of the Gallagher model’s analytic sensitivities, a com- 
parison between those found from the best CST model (adrib = 5) and those from the LST 
model is shown in Table 8.6. Since the maximum deflections differ by 6.7%, it can be 
assumed that all sensitivities would yield closer results if each wing model’s deflection 
behavior were more similar. 

A parametric study to assess the usefulness of analytic sensitivities for future optimi- 
zation usage is performed using the Gallagher model 1 wing under a 100 lb. trailing edge 
tip load. Shape variable xpR is incrementally perturbed to alter the wing planform. The 
trailing edge tip vertical displacement, the second natural frequency, the trailing edge root 
cap stress and the spanwise plane stress CTyy for a centrally located CST wing skin element 
are plotted versus xpR in Figure 8. 1 . First order Taylor series representations for each out- 


put are obtained from 



103 


Tabic 8. 1 - Analytic vs. finite difference xpR sensitivities - Gallagher CST model 1 


Shape design variable: leading edge wing tip x-location (xpj^) 


Output 

Analytic 

Finite Difference Sensitivity 

Parameter 

Sensitivity 

design variable perturbation 
.001 chord .01 chord .1 chord 

trailing edge tip 
displacement 
(in. / in.) 

0.0206 

0.0206 0.0207 0.0213 

1st natural 
frequency 
(Hz. / in.) 

0.991 

0.991 0.994 1.023 

leading edge 
root cap stress 
(psi / in.) 

-269.69 

-269.33 -269.87 -271.61 

leading edge root 
wingskin sigma Y 
(psi / in.) 

-697.16 

-697.33 -697.13 -696.03 




104 


Table 8.2 - Analytic vs. finite difference yR sensitivities - Gallagher CST model I 


Shape design variable: wing tip y-location (yR) 


Output 

Analytic 

Finite Difference Sensitivity 

Parameter 

Sensitivity 

design variable perturbation 
.001 span .01 span .lspan 

trailing edge tip 
displacement 
(in. / in.) 

0.1082 

0.1083 0.1091 0.1173 

1st natural 
frequency 
(Hz. /in.) 

-3.68 

-3.68 -3.63 -3.21 

leading edge 
root cap stress 
(psi / in.) 

537.16 

537.12 535.49 520.00 

leading edge root 
wingskin sigma Y 
(psi / in.) 

590.29 

589.99 588.44 571.43 
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Table 8.3 - Analytic vs. finite difference xjtr sensitivities - Denke CST model 


Shape design variable: leading edge wing tip x-location (xpr) 


Output 

Analytic 

Finite Difference Sensitivity 

Parameter 

Sensitivity 

design variable perturbation 
.001 chord .01 chord .1 chord 

trailing edge tip 
displacement 
(in. / in.) 

7.798x1 O' 6 

7.830x1 0‘ 6 

8.250x1 O’ 6 

9.08 lxl O' 6 

1st natural 
frequency 
(Hz. / in.) 

1.617xl0" 3 

1.567xl0* 3 

1.129x10 3 

0.253x1 0* 3 

leading edge 
root cap stress 
(psi / in.) 

0.180 

0.180 

0.176 

0.139 

leading edge root 
wingskin sigma Y 
(psi / in.) 

0.667 

0.667 

0.660 

0.616 
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Table 8.4 - Analytic vs. finite difference y R sensitivities - Denke CST model 


Shape design variable: wing tip y-location (y R ) 


Output 

Analytic 

Finite Difference Sensitivity 

Parameter 

Sensitivity 

design variable perturbation 
.001 span .01 span 

lspan 

trailing edge tip 
displacement 
(in. / in.) 

1.082x1 O' 3 

1.083xl0' 3 1.089x1 0‘ 3 1.155xl0' 3 

1st natural 
frequency 
(Hz. / in.) 

-1.61 

-1.61 -1.60 

-1.53 

leading edge 
root cap stress 
(psi / in.) 

8.32 

8.32 8.33 

8.40 

leading edge root 
wingskin sigma Y 
(psi / in.) 

5.63 

5.64 5.65 

5.85 
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Tabic 8.5 - Analytic vs. finite difference Aj sensitivities - Gallagher CST model 1 


Shape design variable: leading edge root cap area (Aj) 


Output 

Analytic 

Finite Difference Sensitivity 

Parameter 

Sensitivity 

design variable perturbation 



.001 Aj .01 Aj 

1A, 

trailing edge tip 
displacement 

-0.0716 

-0.0716 -0.0716 

-0.0709 

(in. / in. 2 ) 




1st natural 




frequency 
(Hz. / in. ) 

2.65 

2.65 2.65 

2.63 

leading edge 
root cap stress 

19152.7 

19150.7 19134.2 

18970.2 

(psi / in. 2 ) 




leading edge root 
wingskin sigma Y 

21046.9 

21044.8 21026.6 

20846.4 

(psi / in. 2 ) 







Table 8.6 - CST vs. LST analytic shape sensitivities - Gallagher CST model 1 


Output parameter: trailing edge wing tip z-displacement 
CST nodal displacement = 1 .253 in. 

LST nodal displacement = 1 .343 in. (6.7 % difference) 


Design variable 

Analytic Sensitivity 

Percent 

difference 

CST 

model 

LST 

model 

X FL 

0.044846 

0.049825 

9.99 

*AL 

-0.057923 

-0.064331 

9.96 

*FR 

0.022222 

0.025040 

11.25 

X AR 

-0.009145 

-0.010533 

13.18 

yr 

-0.115513 

-.127102 

9.12 

yR 

0.115513 

0.127102 

9.12 

a 

0.392320 

0.399118 

1.70 



vailing edge root cap avesa (p*i) vailing edge tip displacement (in ) 







( 8 - 2 ) 


/(■W m fl*FR\ 0 ) + 


d x FR 


( X FR 

0 



and also plotted. Here, x FR \ is the original value of xpR and /(• x F r\ q ) is the value of 
any parameter. Additionally, reciprocal first order Taylor series approximations (Ref. I) 
are calculated from 


f( x Fp) ~f( x FR\J (■ r F/?l 0 + a ) 


2 d/ (•*>/?) 
dx FR 


f 1 - 

1 1 

ol X FR + a 

X FR^ 0 + a ) 


(8-3) 


Here, ‘a’ is an offset variable to allow for our = 0 case. Figure 8. 1 shows the 

reciprocal approximation when a=-30. As can be seen, first order approximations to the 
non-linear data yield good accuracy for relatively large perturbations in xp^. 


8.3 Computation time assessment 


The Gallagher model 1 wing is used for evaluation of CPU time required for analytic 
sensitivity calculation. A CPU breakdown of each section of the finite element program is 
shown in Table 8.7 with an explanation as follows. Static solution time includes solving 
for every degree of freedom’s displacements and all finite element stresses. Dynamic solu- 
tion time includes computing all natural frequencies and mode shapes (equal to the num- 
ber of degrees of freedom). Design variable sensitivity time includes calculating all 
displacement, element stress and natural frequency sensitivities with respect to any single 
shape or size type variable. 

In looking at the model with four divisions per section, it can be seen that the total 
CPU time to compute the model’s displacements, stresses, natural frequencies and mode 
shapes is 271.81 1 seconds, with either an additional 14.503 seconds to calculate one set of 
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Table 8.7 - Finite element code CPU breakdown - Gallagher CST model I 


CPU seconds / module 


Program 

Number of dummy ribs per section 

module 

0 

1 

2 

3 

4 

degrees of freedom 

90 

180 

270 

360 

450 

Finite elements 

156 

264 

420 

552 

684 

program initialization 

0.100 

0.152 

0.227 

0.316 

0.363 

form global stiffness 
and mass matrices 

0.098 

0.176 

0.316 

0.426 

0.496 

static solution: 

* displacements 

* stresses 

0.113 

0.016 

0.855 

0.031 

2.656 

0.051 

6.020 

0.062 

11.402 

0.074 

dynamic solution: 

* natural frequencies 
and mode shapes 

3.527 

20.065 

60.336 

138.784 

259.476 

shape variable sensitivity: 

* w.r.L one variable 

* w.r.L all shape variables 

2.695 

18.865 

4.883 

34.181 

8.277 

57.939 

11.292 

79.044 

14.503 

101.521 

sizing variable sensitivity: 

* w.r.L one variable 

* w.r.t. all size variables 

0.046 

7.222 

0.098 

25.846 

0.177 

74.144 

0.263 

145.385 

0.361 

246.719 

solution time: 

* no sensitivities 

3.854 

21.279 

63.586 

145.608 

271.811 

solution time: 

* all sensitivities 

29.941 

81.306 

195.669 

370.037 

620.051 
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analytic shape sensitivities or an additional 0.361 seconds to calculate one set of analytic 
size sensitivities, for a worse case run time of 286.3 14 seconds. Using finite differences, 
this same model would have to be analyzed twice (543.622 seconds total) before even pro- 
ceeding with the differencing calculations, thus showing the huge computational advan- 
tage of computing the sensitivities analytically within the program. 

Notice the disproportionate amount of time required to calculate the complete set of 
model natural frequencies and mode shapes. In the future a new eigenproblem solver will 
be added that will solve for only a user specified number of frequencies which will drasti- 
cally cut down the run time. 



CHAPTER 9 


CONCLUSION 


A fresh examination of wing finite element modeling practices shows that accurate 
displacements and natural frequencies can be obtained using simple triangular elements 
(such as the CST and LST) together with rod elements. Smoothing and averaging of 
resulting stresses lead to globally reliable stress predictors. With automatic mesh genera- 
tion and dummy elements, finite element models of wings, including their skins, ribs and 
spars, can be generated efficiently. The elements used make it possible to obtain deriva- 
tives of behavior functions such as displacement, stress and natural frequency analytically 
with respect to shape and sizing design variables. 

Extensive numerical tests comparing predictors of the current capability developed 
with experiments and commercial finite element codes are described. Analytic sensitivity 
calculations are compared to finite difference results and optimization package usage of 
these sensitivities is explained. Thus, the optimization of wing structural systems during 
conceptual or preliminary design phases can be made practical and computationally cost 
effective. 

Future extensions of this work include: 

a) composite material capability, 

b) efficient computation of low frequency modes, 

c) skin buckling predictions, 

d) integration with aerodynamic loads, 

e) reliable weight estimation for as-built wings. 

The work can also be extended to the modeling of whole airplanes. 
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APPENDIX A 


ELEMENT STIFFNESS, STRESS AND MASS MATRICES 


A.l Rod element 


A. 1.1 Stiffness matrix 


The stiffness matrix for a linear, three dimensional rod in its local coordinates (Ref. 
14) is given by 




(A-l) 


where A is the cross-sectional area, E is the Modulus of Elasticity, L is the element length 
and the two degrees of freedom are the axial displacements uj and U 2 only. To transform 
this to the global system, the equation [kglobal] = [T] t x [klocal] x [T] is used with 


m 


cx cy cz 0 0 0 
0 0 0 cx cy cz_ 


where cx = (x 2' x lV L , cy = ^ y 2' y l V L , cz = (z 2' z iVl (directional cosines) and 


(A-2) 


L = J(x 2 -x,) 2 + (> 2 ->i) 2+ (*2 -z i) 2 


(A-3) 


to arrive at the symmetric 6x6 global stiffness matrix 
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A.1.2 Stress matrix 


The axial stress is a scalar and is found through Hooke’s stress/strain law a = Ee . 
To find the local strain in the rod, this is simply the change in length divided by the origi- 
nal length in matrix form as 

e iocal = = ll - 1 l] { U L> (A_4: 

with {U l } t = {u,,u 2 }- Global strain is found using the previous transformation 
[T]{ U G } in place of {U L }: 


We[-i 0 in {v a } 


(A-5) 


where {Uq} 1 = {uj, vj, wj, u 2 , v 2 , w 2 ). In matrix form the stress can now be given as a 
scalar using known global displacements as 
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(A-6) 


or in explicit form as 


a = 


'[(X 2 -Jt,) (m 2 -«,) + ( v 2 — ) ( v 2~ l *i) + (<2 _;: l) ( w 2“ w |)] 

(x 2 -x l ) 2 + (y 2 -y l ) 2 + U 2 “£|) 2 


(A-7) 


A. 123 Mass matrix 


The mass of a rod element is equal to p AL where p is the mass density with A and L 
defined previously. To form the 6x6 lumped mass matrix in the global system, the mass is 
allocated evenly to each degree of freedom by dividing by the number of nodes. Thus 


m roJ 


100000 
0 1 00 00 
P AL 0 0 1 0 0 0 
2 000100 
000010 
000001 


(A-8) 


A.2 CST element 


A.2.1 Stiffness matrix 


The derivation of a constant strain or constant stress triangular element is taken 
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directly from Ref. 23 and Figure A. I. Its basic assumptions are: 

1 . isotropic material 

2. uniform thickness ‘t’ 

3. plane stress state 

4. constant strain in field 

Based on element geometry, the displacement state in local coordinates is 


« (*,y) 


(- (b - s) x-hy) u p + (-r(x-r) +h (y -r) ) u Q + xbu R 

y (x, y)_ 


(-( b-s)x-hy)v p + (-s(x-s) +h(y-s))vQ + xbv R 


where b, s, h and a are local geometric variables (b is the major base, s is the minor base, h 
is the element height and a is the total area). The strain-displacement relation obtained by 
differentiating the above with respect to x, y and z is then 


e xx 1 

II 

SI- 

-(b-s) 

0 

-s 0 


0 

-h 

0 h 

e * 

*y 


-h 

1 

1 

h -s 


0 0 


= [s]{^} 


(A-10) 


and the displacement transformation law from global coordinates is 
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[A] 1 U<;) 


where 

/i m, n, 
l 2 m 2 n 2 



and 



000 
0 0 0. 


(A- 1 1 ) 


(A- 12) 


(A- 13) 


The 1, m and n terms are direction cosines of the local axes with respect to the global 
axes with 


1 s 

l \ = r(*/r(r) (x Q -*p) -Xp) 


(A- 14) 
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1 s 

= {y Q~ y r ) ~ y P ) 

(A-15) 

= ft ( z r~ (ft} ~ z r^ ~ z p) 

(A- 16) 

II 

1 

* 

(A- 17) 

= £< v cr»> 

(A- 18) 

i , , 

- g Uy Z P ) 

(A- 19) 


Due to the plane stress assumption, Hooke’s Law gives 
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(A-20) 


where E is the Modulus of Elasticity and u is poisson’s ratio. To find the stiffness matrix 
in local coordinates, integrating over the area via 


M = m T im*v = #]"&>][*]<*• - [i,]+[*,] (A-2i) 

V S 


V. 


assuming a constant thickness ‘t’ results in 
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and 
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(A-23) 


The 9x9 global stiffness matrix for the CST is then 


M - M MM 


(A-24) 


A.2.2 Stress matrix 


Using Hooke’s Law again, the 3x 1 CST stress vector is given by 



= [d] [b] { U L ] = [d][b][a]{U g } 


(A-25) 
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A. 2.3 Mass matrix 


The mass of a CST element is equal to pAr where p is the mass density with A and t 
defined previously. To form the 9x9 lumped mass matrix in the global system, the mass is 
allocated evenly to each degree of freedom by dividing by the number of nodes. Thus 


100000 
0 10 0 0 0 
_ 0 0 1 0 00 
3 000100 

000010 
000001 


(A-26) 


A.3 LST element 

A.3.1 Stiffness matrix 

The derivation of a linear strain triangular element is taken directly from Ref. 17 and 
Figure A.l. It’s basic assumptions are: 

1 . isotropic material 

2. uniform thickness 

3. plane stress state 

4. linear strain in field 

The local 12x12 stiffness matrix is given by 


= H r HH 


(A-27) 
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(A-28) 


(A-29) 


where c n = c 12 = % . v 2 )> c 12 = vE/ (l- v 2 ) ^d c 33 = E ° * V)/ 2(l - v 2 ) = E/ 2(l + v)- 7116 
global geometry variables { B } = {aj, clj, bj, b 2 , b 3 ) are linked to the local geometry 
variables {G} = { b, s, h, a} by 
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a | = -h 32 — h aj = 0 
b | = s - b b2 = -s bj = b 


To find the IXx IX global stiffness matrix, transformation is the same as for the CST in 
section A. 2. 1 using 


r. i _ 
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(A-30) 


A.3.2 Stress matrix 


Using Hooke’s Law, the 9x1 stress vector consisting of <J XX . Gyy and (J xy for each of 
the three end nodes P, Q and R is 
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(A-31) 
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(A-34) 
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CST clement [P Q R] 


Yjocal 



P 


LST element [P Q R T U V] 


Figure A. 1 - Triangular membrane elements used 


APPENDIX B 


ELEMENT COORDINATE SHAPE VARIABLE DERIVATIVES 


B.l d{X}/d(x FL , x AL , x FK , x AR , yp, y R ) 

Based on the geometry in Figure 3. 1 . it can be seen that every ‘y* value at point 
T is a linear combination of yp and y R such that 


v, = y L + P s (>/ ? - y L ) = (>-/»,) >’l + P*yR 


(B- 1 ) 


where p s = percent span ratio in the y-direction and is given by 


Ps 


- v < ~ - v z. 

y R~y l 


Now, differentiating yj with respect to the six shape variables yields 


(B-2) 


dy/dxpL = 0 
dy/dx^ = 0 
dy/dxpR = 0 
d ys /dx AR = 0 

dy/dyL = 1 - Ps 
dy/dy R = p s 


For the V values at point ‘i,’ if ‘i’ is along either the wing root or wing tip, the situa- 
tion is the same as for the *y* values above. Along the root, V is given by 
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= X FL + Prc ( X AL “ X FL> = < 1 “ /»«■) X FL + Prc X AL 


(B-3) 


while along the tip, V is given by 


x i = x FR + Pn ( x AR~ x F^ = 0 ~P,^ x FR + Pu x AR 


(B-4) 


where p rc = percent chord ratio along the root and p, c = percent chord ratio along the tip 
and are given by 



and 


(B-5) 



(B-6) 


Therefore, differentiating Xj along the root with respect to each shape variable gives 


dXj/dxpL = 1 - p rc 
dx/dx^ = p rc 
dx/dxp-R = 0 
dxj/dx^ = 0 
dx/dyL = 0 

dxj/dyR = 0 


while doing the same along the wing tip yields 
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dx/dxpL = 0 

dx/dx^L = () 

dx/dxpR = I - p tc 
dx/dx^ = p, c 


dx/dy L = 0 

dxj/dy^ = 0 


For all remaining nodes, the ‘x’ values are a linear combination of all four ‘x’ -valued 
shape variables. Fortunately, due to the nature of the wing geometry, since we know the 
‘x’ derivatives along both the root and tip, it is a straightforward process to interpolate 
what they should be for any point ‘i’ across the span. In other words 

X: = (x) ,+ Pc [(*;), - (*.-) ,] (B-7) 

1 1 root ' s t tip 1 root 

or 

*, = 1(1 -P r J X FL + Prc X A J O “/»,) + [(1 ~ P, c ) X FR + ( B ' 8 ) 

where p s has been defined above (eqn. B-2). Differentiating with respect to each shape 
variable then yields the more general and final form of 


dx/dxpL = (1 - p rc )(l - p s ) 
dx/dx,^ = p rc (l - p s ) 
dx/dxpR = (1 -p,c)p s 
dXj/dxAR = PtcPs 
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dxj/dy L = 0 
dx/dy R = 0 


As one can see, this reduces to the simplified forms along the root and tip above when 
p s equals 0 and 1. respectively. When the depth distribution is given in global coordinates, 
then all ‘z’ values at point ‘i’ are independent of these shape variables so that their deriva- 
tives are equal to zero. If the depth distribution is dependant on the wing trapezoid shape, 
sensitivities with respect to shape variables must be included. 

In summary. 



(1 -P rc ) U -P s ) 
0 
0 


(B-9) 
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(B- 10) 


a f*« ] f (i-PtJp, 

^ll;H s 


(B- 11) 



(B- 12) 
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(B- 13) 



(B-14) 


B.2 d{X}/d(a) 

Based on the geometry in Figure 3. 1, it can also be seen that for any nodal point ‘i\ 
if given a ‘y’ value, then the corresponding ‘x’ value is given by 

x i = y-tana + C (B- 15) 

where C is any constant. If all ‘y’ and ‘z’ coordinates are assumed to be independent of the 
sweep angle, differentiating with respect to the design variable Ot gives 


3 f x i ] f y t ( sec (a)) 2 

5^ \ y i = 0 

1 Zi 1 0 


(B-16) 



APPENDIX C 


SHAPE VARIABLE SENSITIVITIES 


C.l Global displacement sensitivity with respect to any shape variable 


From the basic static equation f K] { U } = { F } . one can differentiate with respect to 
any shape design variable (3 to get 


M 


3{t/} 

ap 


at/n 
+ ap 


{U} 


d{F} 

ap 


(C-l) 


where [K] is the global stiffness matrix, { U } is the global displacement vector and { F } is 
the global load vector. 

For any loading case in which the applied loads are independent of model geometry, 


a{t/} 

ap 


= -m 


_,am 

~W 


{£/} 


(C-2) 


With [K] and {U} having already been computed, and the partial derivative of 
displacement with respect to any shape design variable desired, it is only necessary to 
compute the global stiffness matrix derivative. This is done on an element by element 
basis and the individual results are then merged as done when forming [K] previously. 

C.1.1 Rod element stiffness sensitivity 


Using chain rule differentiation, the derivative of a rod element stiffness matrix with 


respect to any shape design variable is 
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3 [*<ll _ 8 C*<il 3{->f} , rr 

“5F a{X)~5F ' 

where {X } = { Xj , yj, Z|, ^ Y2- ’ L i\- The partial differentiation of {X} with respect to any 
design variable has previously been calculated in Appendix B. To find the partial deriva- 
tive iif [k^] with respect to the rod element’s nodal coordinates { X J. straight-forward 
chain rule differentiation is carried out (with the following simplifications: Ax = X 2 - Xj, 
Ay = y 2 -yj, Az = Z 2 - Z|) with ‘i’ ranging from 1 to 6: 



~[AA] i ~[AA)] 

ax, 

~[AA) i [AA) l 


(C-4) 


where 
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Then, the 6x6 rod element global stiffness matrix sensitivity with respect to shape design 
variable p is 

3 [*e] 3 [*«] Ml . 3 [*c] 3AT 2i fc-II 

ap " axj ap ax 2 ap ax 6 ap 


C.i.2 CST element stiffness sensitivity 


Chain-rule differentiation of 9x9 CST stiffness matrix [k^] with respect to any shape 
design variable gives 


a [*c] _ a [*c] d{X} 

ap 3 { xj ap 


(C-12) 



I3S 


so that only the partial derivative of [kg] with respect to {X} needs to be found. To calcu- 
late this, differentiation of the matrix expression for [k<j] yields 


1 -1 




^ = ' A ' r '‘Jn3o MA|, aw |A| ^w 


I*, I IA1 


(C-13) 


where { X } = { x j , yj, zj, x 2 , y 2 , z 2 , x 3 < Y 3 > z 3 ^- All undifferentiated matrices are known so 
that the only unknowns are the transformation matrix derivatives and the local stiffness 
matrix derivatives each with respect to nodal coordinates. 

Before proceeding, all geometric variables will be linked to each other through 
Figure A. 1 and the following equations: 

{L} = { lj, 1 2 , 13 } = function of {X} only where 
1 , = [(x 3 - x 2 ) 2 + (y 3 - y 2 ) 2 + (z 3 - z 2 ) 2 ] 1/2 

1 2 = [(x 3 - xj) 2 + (y 3 - y ^ 2 + (z 3 - zj ) 2 ] 172 

1 3 = [(x 2 - xj ) 2 + (y 2 - yj ) 2 + (z 2 - zj ) 2 ] 1/2 
{G} = {b, s, h, a} = function of { L} only where 

b = l 3 

s = (1 2 2 + 1 3 2 - 1 j 2 )/(21 3 ) 
h = [l 2 2 - d2 2 + 13 2 - 1i 2 ) 2 /(41 3 2 )] 1/2 
a = ( J / 2 ) 1 3 [1 2 2 - (1 2 2 + 1 3 2 - 1j 2 )/(41 3 2 )] 1/2 
[kj = function of Young’s Modulus, thickness and {G} only 
[A] = function of {G} and {X} only 
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C. 1.2.1 d[k L ]/d{X}: 


Chain rule differentiation of the 6x6 [klocal] with respect to vector {X} gives 

a [*t] _ %] 3(C( 3{L) 

3{X} d{G}JJTTd{X} 


where 
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(C-2 1 ) 
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(C-22) 
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along with 









a{G} 

d{L] 


0 

-1 

h 

2hl\ 
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8a 


0 

h 

h 

/ 2 (/?+/H 2 2 ) 

2/1/3 

i 2 (i\ + i\-i\) 
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(C-24) 


and 
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d{L} 

aw 


0 0 0 


x 2 *3 ^2 - V 3 z 2 <-3 x 3 X 2 - v 3 - V 2 ^3 Z 2 


x l x y .' i ^3 H ^ 3 


/, 


2 ^2 ^2 ‘2 
x \~ x 2 X \~ X 2 Z \~ Z 2 X 2 ~ X l >’ 2 “- V l Z 2 ~ Z \ 


/, /, /, 


0 0 0 




'1 


■*3 r | - v 3 V 1 ‘•3 2 1 


/, /, ly 1 3 I) I ) 


1 2 
0 


l 2 l 2 

0 0 


(C-25) 


Thus, to find the derivative of [k L ] with respect to any X v chain rule summation yields 


_ d [ k i ac, aw ac 4 a {/.} 

^ ” ^alw ~a^ + '' + aG 4 a w ax, 


(C-26) 


where 



a{L> 


is the 1x3 row *j* of 


a{G> 

aw 


and 


aw 

a*, 


is the 3x1 column T of 


aw 
aw ' 


d[A]/d[X]: 


Chain rule differentiation of the 6x9 transformation matrix [A] with respect to {X} 
gives 

»[A] _ a W a{G} aw a[A] 
d{x } a{Gjawaw aw 


(C-27) 
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d[a] 

where ^ is the total derivative since there is explicit dependance of [ A] 
The following new derivatives are used: 


a[A] 

3G; 


[AA], [0] 10] 

[0] [AA) i [0] 
[0] [0] [AA] J 


(i = 1 to 4) 


with 


[AA ] . = 


[AA] 0 = 


[AA], = 


5 (X 2 -X,) 

•s ( V 2 _ >1 ) 

^ (^2 ) 

b 2 h 

b 2 h 

b 2 h 

(x 2 -x 1 ) 

(yi~y i) 

(z 2 ~*i) 

b 2 

b 2 

b 2 -J 

(x 2 -x l ) 

(yi-y i) 

(z 2 -Zi) 

bh 

bh 

bh 

0 

0 

0 


-JTj) -x, 

*3 -fa 

h 

2 

h‘ 

0 


0 


h 1 

0 


[AA] 4 = 


000 

000 


on {X}. 


(C-28) 


(C-29) 


(C-30) 


-Z| 


(C-31) 


(C-32) 
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and 



[BB] i [0] [0] 

[0] |5S], [0] 

[0] [0] [BS],. 


(i = 1 to 9): 
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\_ 

b 


0 0 


(C-33) 


(C-34) 
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(C-35) 


[BB ] 3 = 
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r 00 

b 


(C-37) 
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bh 
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0 

0 
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s 

bh 

[ 

b 


(C-39) 


[BB] y 



0 0 OJ 


(C-40) 


[BB} % 


0 y 0 

h 

0 0 0 


(C-4 1 ) 


[BB) 9 = 


o i o 

h 

0 0 0 


(C-42) 


C.13 LST element stiffness sensitivity 




Chain-rule differentiation of [k^] with respect any shape design variable gives 


a[* c ] _ a [*c] 3{x> 
lijr “ dTx] ap 


(C-43) 


so that only the partial derivative of [k^] with respect to {X} needs to be found. To calcu- 
late this, differentiation of the LST matrix expression for [kc] yields 
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apU - r ^ a . r 9 W - 

Tm = [A| ly 5W t|A| 5W |A|+ 5W '*4 ,A| 


where {X} = {xj, yj, zj, x 2 . y 2 , z 2 . x 3 . V 3 } even though the LSThas twice the number of 
nodes of the CST element. The reasoning is that since the side nodes are assumed to be 
placed at the mid-point of each side, their location depends on the comer nodes. All undif- 
ferentiated matrices are known so that the only unknowns are the transformation matrix 
derivatives and the local stiffness matrix derivatives with respect to nodal coordinates. 

Like the CST element, all geometric variables will be linked to each other through 
Figure A. 1 and the following equations and vectors: 


{ L } = { 1 1 1 2 I 3 } = function of { X } only where 
lj = [(*3 - * 2 ) 2 + (Y3 ’ y2) 2 + ( z 3 ‘ z 2) 2 ] 1/2 

i 2 = t( x 3 * x i) 2 + (y3 ■ yi) 2 + ( z 3 ■ z i) 2 ] 1/2 

1 3 = [(x 2 - x,) 2 + (y 2 - yj) 2 + (z 2 - zj) 2 ] 1/2 

{G} = {b s h a} = function of {L} only where 
b = 1 3 

s = (1 2 2 + 1 3 2 - 1, 2 )/(21 3 ) 
h = [I 2 2 - d 2 2 + b 2 ~ 1 i 2 ) 2 /(41 3 2 )] 1/2 
a = (V 2 ) 1 3 [1 2 2 - (1 2 2 + 13 2 - 1i 2 )/(41 3 2 )] 1/2 
(B} = {al a2 a3 bl b2 b3} = function of {G} only where 
al = -h 
a2 = h 
a3 = 0 
bl = s - b 
b2 = -s 
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b3 = b 

[kj = function of Young’s Modulus, thickness, {G} and { B } only 
f A] = function of { G } and { X } only 

d[k L ]/d{X>: 


Chain rule differentiation of 12x12 [kjJ with respect to {X} gives 


a {B} d{G} d[k L ] 9{G) 

aixT " + 


From the CST element in Appendix C.1.2 we have the 4x9 derivative matrix d ^/<j{X) 
already. Differentiating each component of { B } with respect to each component of {G } 
gives the 6x4 matrix 


0 0-10 
0 0 10 
d W 0 0 0 0 

d{G) -110 0 

0-100 
10 0 0 


(C-46) 


All that remains is the single derivative of [kj with respect to both {B} and {G} 
where [kj = [M] T [N][M] from Appendix A.3.1. Differentiation of the 12x12 local stiff- 
ness matrix against the six components of vector {B} yields 


3{B> 


m T m 


3 [a/] 


3[w] 

3{B) 


T 

[H] [M] 


{CA1) 


with all undifFerentiated matrices previously known ([N] is not a function of { B }). 



I4K 


Straightforward differentiation of 9x12 [M] with respect to { B } yields 




I 

2 a 
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3 0000000000 
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d[Af] 
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0 

00004000 
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J 0000 0 -10000 04 
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0000-1 0000040 
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(CMS) 


(C-49) 


(C-50) 
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0 0000000000 
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( C - 5 1 ) 
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2 a 
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( C - 53 ) 
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Differentiation of the 12x12 local stiffness matrix against the four components of 


vector { G } gives 


d [ k i 

d{G} 


[A/] r [ /V] 


a&f] 
ai c} 


+ 


a [a/] 
ajc7 


[N] [A/] 


(C-54) 


with all undifferentiated matrices previously known t[N] is not a function of {G}). 
Straightforward differentiation of 9x 1 2 [M] with respect to {G } is simplified since 
[M] is not a function of Gj, Gi or G 3 . Therefore, 


a [a/] _ a [a/] _ a[A/] 
5G, dG 2 aGj 

and 


(C-55) 


a [a/] 

aG7 


3*. 
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-h 
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-b 3 
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4b 2 
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4b 3 
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3b 2 

0 
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0 

Ab x 
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4*3 

0 

0 

0 


0 
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0 

3^3 
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0 

0 

4£» 2 
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4 *, 

0 

0 

3 fl i 

0 

~ a 2 

0 

-a 3 

0 

4a 2 

0 

0 

0 

4a 3 

0 


0 

3 a 2 

0 

" a 3 

0 

4a, 

0 

4a 3 

0 

0 

0 

~ a i 

0 

-a 2 

0 

3a 3 

0 

0 

0 

4a 2 

0 

4a, 

3<*1 

3 b x 

-°2 

-b 2 

-a 3 

"*3 

4a 2 

4b 2 

0 

0 

4a 3 

4b 3 



3a 2 

3b 2 

-a 3 

— *3 

4a, 

4 *, 

4a 3 

4£> 3 

0 

0 




~b 2 

3a 3 

3b 3 

0 

0 

4a 2 

Ab 2 

4a, 

4*,. 


(C-56) 


Thus, all necessary derivative matrices are known and ^lV^x] can calculated. 
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d[A]/d[X]: 

Since the 12x18 LST transformation matrix is formed from two 6x9 CST transforma- 
tion matrices, the derivative with respect to {X } is simply the CST’s transformation deriv- 
ative with respect to {X} (Appendix C. 1.2) used twice as 


d A 

JJX} 


a[A] 

9{X} 

[ 0 ] 


[ 0 ] 

a[A] 


a{xj 


(C-57) 


C.2 Global stress sensitivity with respect to any shape variable 

Unlike displacement sensitivities, stress sensitivities can be calculated on an element 
by element basis. 

C.2.1 Rod element stress sensitivity 

The derivative of the rod element’s scalar axial stress is 

Bo Bo 3{X} Bo B{U C } lnm 

5P"S{X) + 3{(/ c } 5T" 

where (X} = {xi, y j, Zj, x 2 , y 2 , z 2 } and Uq = {ui, vj, w^ u 2 , v 2 , w 2 }.With the design 
variable displacement and coordinate derivatives having been previously calculated, all 
that is necessary is the derivative of the stress equation with respect to {X} and {Uq}. For 
completeness, this results in 
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= — | l 2 («-, - «,) - 2 [ (« 2 -u,) + (y : -y,) ( v 2 ~ v i) + ( z 2~ z i) (*v 2 - w, ) i (x,-x,)| (C-59) 


L 4 


^ \L 2 (V 2 - V,) - 2 [ (Jtj-Jt,) (« 2 ~«,) + (V : -V,) + (Z-. ~ Z, ) (*V 2 -«,)| (>;->•,)) (C*60) 

dy i l a 


^ p 

= _ ( i. 2 (« , - vv, ) - 2 I («-,-«,) + (>•;->•,) < v : ~ *’| 

3; i L 4 


:,) (C-61) 


(C-62) 


(C-63) 


(C-64) 


do E 


= ~~2 ( x 2 ~ x l '> 


d“i L 


(C-65) 


do E . , 

5 


(C-66) 


E < ^ 


(C- 67) 


(C-68) 


(C-69) 


(C-70) 
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C.2.2 CST element stress sensitivity 


Using chain-rule differentiation on the CST stress equation yields 


cl 

dP 


o 

.1. 

o 

y 


= [«] H + &] { u k ] + H { t/ « > 




(C-7 1 ) 


where {U G } = {uj, vj, w lt u 2 , v 2 , w 2 , U3, V3, w 3 } and all undifferentiated matrices are 
previously known. The displacement sensitivity vector is also known as it was calculated 
above. Therefore, to find the CST stress sensitivities with respect to shape, only the trans- 
formation matrix derivative and [B] matrix derivative are needed. Fortunately, the trans- 
formation derivative has already been found to be 


a[A] 

~W = 


a W d{G} 8{L} d [A] 1 8{X} 

8 {G} 8 {L} a {X} 8 {X} 


(C-72) 


Thus, to find d ^ B V d (p), use the chain rule to get 


a [#] _ a [ b] a {G} dm a{x> 

“ d{G}J{lJJ{X}~W' 


Here, the only unknown is dfB V d { G } which can be explicitly found as 
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B[B) 
dG j 



0 00 


X 0 -X 0 0 


~~ 4- 00 

b 2 h b 2 b 2 h 


(C-74) 


BIB) 

BG2 


u 0 -u 0 00 

0 0 0 0 0 0 

,° Vk 0 -u 00 _ 


(C-75) 


3[fi] 

ac 3 



o 

o 

(b-s) 

bh 2 


bh 2 
0 

0 


0 

0 

s 

bh 2 




(C-76) 


a [ 5 ] 0 0 0 o o 0 

-gQ- = 0 0 0 0 0 0 

4 [0 0 0 0 0 0 _ 


(C-77) 


C.2.3 LST element stress sensitivity 

Using chain-rule differentiation on the LST stress equation yields the 9x 1 stress 
derivative vector with respect to any shape design variable (5 as 
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(C-78) 


where all undifferentiated matrices are previously known. The displacement sensitivity 
vector is also known as it was calculated in Appendix C. 1. Therefore, to find the LST 
stress sensitivities with respect to shape, only the transformation matrix derivative and 
[M] matrix derivative are needed. Fortunately, the transformation derivative matrix has 
already been found above to be 



a[A] 

d{X} 

[ 0 ] 


[ 0 ] 

a[A] 


aw 


aw 

~W~ 


(C-79) 


Thus, to find d[M]/d((5), use the chain rule as before to get 

a [M] _ ,d[M]d{B } a{G} aw a[M] aw aw aw 
" ( awa{Gj aw aw + a{G} aw aw ,_ 5F 


All of these derivative matrices have been previously calculated in Appendix C. 
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C.3 Natural frequency sensitivity with respect to any shape variable 


Differentiating the eigenvalue equation 
[K-v> 2 M] {<D} = {0} 


with respect to any shape design variable P yields 


9 ( 0 ? 

w = 


♦f 


d[K] ,d[M] 


L dP 


COT 


dP 


♦. 


\M] <t>; 


(C-8 1) 


(C-82) 


which is actually the sensitivity of the eigenvalue X = to 2 . The global mass sensitivity 
matrix ^ is the only new entry. To calculate this, individual element mass sensitiv- 

dP 

ity matrices are calculated and then merged similar to what is done for the global stiffness 
matrix. 

C.3.1 Rod element mass matrix sensitivity 


Differentiating the 6x6 global rod mass matrix (Appendix A. 1.3) yields 


d [M rod \ _ d[Wja{X} 
5]j 9{X} 


(C-83) 


0 r x] 

where — has been derived previously and {X} = {xj, yj, zj, X 2 , y 2 , Explicit dif- 
ferentiation of with respect to {X} yields (since it’s length L is a function of {X}) 
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P‘4 

ax t 

2L 


pA 

dX 2 

2L 

a («,J 

P A 

ax 3 

2L 


= £d, 

ax 4 

2L 

3I*„J 

- 

3x 5 

2L 

3t«rJ 

= 

3X 6 

2L 


-2 <- 1 - 




(C-84) 


(C-X5) 


(C-86) 


(C-87) 


(C-88) 


(C-89) 


where p is the density, A is the cross-sectional area, L is the length and I is a 6x6 identity 
matrix. 


C.3.2 CST element mass matrix sensitivity 


Differentiating the 9x9 global CST mass matrix (Appendix A.2.3) yields 
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d[A/^y] 3 3|X} 


w 


a {x} dp 


(C-90) 


where 


d{*} 

dP 


has been derived previously and {X} = {xj, yj, Z|, X 2 , y 2 . z 2 . x 3 - >’ 3 < z 3 l- 


Chain rule differentiation of Mqjt w ‘ t * 1 res P ect to (X) yields (since it’s area ‘a’ is entry 
#4 in the geometry vector { G ) = { b, s, h. a } ) 


d d d {G} d{Z/} 

d{X} = d{G} dUTdTxT 


where {G } = {b, s, h, a} and {L} = { 1 j T 1 2 , 13 } . Since M CST is not a function of b. s or h, a 
more exact form gives 


d \.M C st\ d[M^_jy] da d{L} 
d{X} = da dTHdW 

d{L} 

with 3x9 matrix ^ already having been calculated. Thus 


d[M CS j] _ pr dG 4 d{L} 
— dX^ TdT ~L)~dX~ [ 


(C-93) 


where p is the density, t is the cross-sectional thickness and I is a 9x9 identity matrix. 



APPENDIX D 


SIZE VARLABLE SENSITIVITIES 


D.l Global displacement sensitivity with respect to any size variable 


From the basic static equation [K]{ U} = { F } . one can differentiate with respect to 
any size design variable K to get 


M 


d{U] , 3[Jf] 


3k 


+ 


9k 


{^1 


d{F) 

9 k 


(D-l) 


where [K] is the global stiffness matrix, { U } is the global displacement vector and { F } is 
the global load vector. 

For any conservative loading case in which the applied loads are independent of 
model geometry. 


d{U} 

~^T 


= -w 


am 




(D-2) 


Again, [K] and {U} are known and the stiffness matrix derivative must be formed via 
merging element by element. But, since a sizing variable exists for each finite element, 
the derivative of the global stiffness matrix of the system with respect to any one size 
variable reduces d ^ K V ( j( 1C ) to a matrix whose only entries are that of the element’s stiffness 
matrix derivative with respect to it’s own size variable. When this extremely sparse matrix 
is multiplied by the corresponding entries in {U}, the global displacement derivative is 
easily calculated. 
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D.1.1 Rod element stiffness matrix derivative with respect to it’s area 


Since [Icq] for rod element ‘i’ is a linear function of it’s area. 




cx 2 

cxcy 
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-cxcy -cxcz 

3IM, 
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- cxcy 
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cx 2 

cxcy cxcz 



-cxcy 

-cy 2 
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cy 2 cycz 



-cxcz 

- cycz 

- cz 2 

cxcz 

cycz cz 2 _ 


where E is the Modulus of Elasticity, L is the element length and cx, cy and cz are the 
direction cosines given in Appendix A. 1 . 1 . 


D.1.2 CST element stiffness matrix derivative with respect to it’s thickness 


Since [k^] for CST element ‘i’ is a linear function of its thickness. 


[k L ] 7-9 [^ + M ir 

Sir = W,-3Z-H,- W, — 57 — M, 


(D-4) 


where 
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(D-5) 


and 
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(D-6) 


k N and k s are the normal and shear stiffness matrices. 


D.1.3 LST element stiffness matrix derivative with respect to it’s thickness 


Since for LST element ‘i’ is a linear function of its thickness. 


a (M , 

Si- 


r-n r3 I*J,' 


JS, 


3t: 


AJ, 


where 


a,, 


= M, sr 



(D-7) 


(D-8) 
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and 
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D.2 Global stress sensitivity with respect to any size variable 

As with the stress shape variable sensitivities, the stress sizing variable sensitivities 
are done on an element by element basis as follows: 

D.2.1 Rod ‘i’ stress sensitivity with respect to rod ‘j’ area 

Differentiating the rod element stress equation with respect to any rod element area 
gives 


= 3(t/ G > 

dAj d{U G } dAj 


(D-10) 


where the displacement derivative with respect to area ‘j’ has been found previously. To 
find the stress derivative with respect to its nodal coordinates, employ straightforward dif- 
ferentiation. This result has been calculated in Appendix C.2.1. 
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D.2.2 Rod ‘i’ stress sensitivity with respect to CST or LST ‘j’ thickness 


Differentiating the rod element stress equation with respect to any CST or LST ele- 
ment thickness gives 


dtj d{U G ) dtj 


(D-ll) 


where the displacement derivative with respect to thickness ‘j’ has been found previously. 
The stress derivative with respect to its nodal coordinates, has been calculated in Appen- 
dix C.2. 1. 


D.2.3 CST ‘i’ stress sensitivity with respect to rod ‘j’ area 


Differentiating the CST element stress equation with respect to any rod element area 
gives 




(D-12) 


where the displacement derivative with respect to area ‘j’ has been found previously. 
Since [D], [B] and [A] are known, the derivative is easily found. 


D.2.4 CST ‘i’ stress sensitivity with respect to CST ‘j* thickness 


Differentiating the CST element stress equation with respect to any CST or LST ele- 
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ment thickness gives 




(D-13) 


where the displacement derivative with respect to thickness ‘j’ has been found previously. 
Since [D], [B] and [A] are previously known, the derivative is easily found. 


D.2.5 LST T stress sensitivity with respect to rod ‘j’ area 


Differentiating the LST element stress equation with respect to any rod element area 
gives 
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where the displacement derivative with respect to area 4 j’ has been found previously. 
Since [C],[M]and [A] are previously known, the derivative is easily found. 

D.2.6 LST ‘i’ stress sensitivity with respect to LST ‘j’ thickness 

Differentiating the LST element stress equation with respect to any CST or LST ele- 
ment thickness gives 


dr 
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d{U G } 

B tj 


(D-15) 


where the displacement derivative with respect to thickness ‘j’ has been found previously. 
Since [C ] , [M] and [A] are previously known, the derivative is easily found. 

D.3 Natural frequency sensitivity with respect to any size variable 


Differentiating the eigenvalue equation 
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[K-w 2 M] {d>} = {0} 


(D-16) 


with respect to any size design variable K yields 


B co 2 

~37 = 




B[K] 

Bk 
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(D- 17) 


which is actually the sensitivity of the eigenvalue A, = co 2 . The global mass sensitivity 

matrix _ is the only new entry. The global stiffness sensitivity matrix has been 
Bk 

detailed in Appendix D. 1 and the global mass sensitivity matrix has the same properties in 
that it’s derivative with respect to any i’th size variable is just the i’th individual mass 
matrix derivative with all other entries equal to zero. 


D.3.1 Rod element mass matrix sensitivity with respect to it’s area 


Since M rod for rod element ‘i’ is a linear function of it’s area, differentiating the 6x6 
global rod mass matrix (Appendix A. 1 .3) yields 


d[A* r< J _ pL 

9A,. 2 


(D-18) 


where p is the density, L is the length and I is a 6x6 identity matrix. 
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D.3.2 CST element mass matrix sensitivity with respect to it’s thickness 


Since M CS t for CST element V is a linear function of it’s thickness, differentiating 
the 9x9 global CST mass matrix (Appendix A.2.3) yields 

d jjjWj _ P^|/j ( D - 1 9 ) 

9r ( 3 1 

where p is the density, A is the area and 1 is a 9x9 identity matrix 


V, . 



